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Summary 



We investigate statistical properties of several classes of periodic billiard models which can be 
regarded as diffusive. We begin by motivating the study of such models in Chap. [T] and reviewing 
how statistical properties arise in Chap.|2l 

In Chap. [3] we consider a periodic Lorentz gas satisfying a geometrical condition, for which 
diffusion has been rigorously proved. We discuss how to estimate diffusion coefficients from 
numerical data and then study their geometry dependence, finding a qualitative change in the shape 
of curves as one parameter is varied. We discuss the application of a random walk approximation 
of the diffusion coefficient and a related Green-Kubo formula. We also consider the effect on the 
diffusion coefficient of reducing the geometrical symmetry. 

In Chap. |4]we study the shape of position and displacement distributions, which converge to 
a normal distribution by the central limit theorem. We find a fine-scale oscillation in the densi- 
ties which prevents them from converging pointwise to Gaussian densities, and relate this to the 
geometry of the billiard domain, giving an analytical expression for the fine-structure function. 
We provide strong evidence that, when demodulated, the densities converge uniformly to Gaus- 
sians, strengthening the standard central limit theorem, and we find an upper bound on the rate of 
this convergence. We further consider the effect of a non-constant distribution of particle speeds, 
showing that the limiting position distributions can be non-Gaussian. 

Chap. |5] treats polygonal billiard channels, where few rigorous results are known. We provide 
numerical evidence that normal diffusion can occur, and that the central limit theorem can be 
satisfied. We develop a picture of how normal diffusion can fail if there are parallel scatterers, 
and we characterise the resulting anomalous diffusion, as well as the crossover from normal to 
anomalous diffusion as such a geometrical configuration is approached. 

In Chap. |6] we extend our methods to a three-dimensional periodic Lorentz gas. We present 
a model with overlapping scatterers exhibiting normal diffusion in a certain regime. Outside this 
regime we provide evidence that the type of holes present in the structure strongly influences the 
statistical properties, and show that normal diffusion may be a possibility even in the presence of 
cylindrical holes. 

We finish in Chap. |7] with conclusions and some directions for future research. 
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Chapter 1 



Introduction 



1.1. Motivation: dynamics of fluids and statistical mechanics 

The dynamics of fluids is comparatively well understoocj^ at the length scales of everyday ex- 
perience: we call this the macroscopic level. At this level, the dynamics is described by partial 
differential equations modelling the behaviour of continua. But we can also think of fluids as made 
up of a vast number of microscopic particles which are orders of magnitude smaller than typical 
macroscopic lengths, and whose interactions we again understand well. The question then arises 
to relate these two levels of description of the same substance. Since many states of the micro- 
scopic system correspond to one state of the macroscopic system (specified by a few macroscopic 
variables such as temperature and pressure), the relation must be of a statistical nature: we seek to 
relate averages over microscopic states to macroscopic phenomena, where there is a separation of 
length scales between microscopic and macroscopic. 

1.1.1. Hard-sphere fluids 

At a microscopic level, the simplest physical^ microscopic picture of a fluid is a collection of a 
very large number of identical atoms or molecules moving through empty space and undergoing 
collisions with each other. Restricting to a classical (as opposed to a quantum-mechanical) descrip- 
tion, we consider a Hamiltonian system with inter-particle interactions described by a short-range 
potential. 

The simplest potential, giving the most naive picture of a fluid, is the hard-sphere potential, 
which jumps from to oo at a distance a; the dynamics then corresponds to spheres of radius a 
in free motion which undergo elastic collisions when they meet. Simulations of a relatively small 
number of hard spheres (hard-sphere fluids) show behaviour which, in certain regimes, resembles 

'The major exception to this is the class of turbulent phenomena. 

^Many seemingly less-physical microscopic models satisfying certain conservation laws give rise to fluid behaviour 
at a macroscopic level. A particularly good example is that of lattice-gas cellular automata, where particles jump 
between neighbouring points of regular lattices and obey certain microscopic conservation laws: see e.g. IRBOll . 
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that of real fluids. The hope is that studying this relatively simple system should give insight into 
the relation between microscopic dynamics and macroscopic fluid dynamics. 

1.1.2. Transport properties 

One of the key aims of (non-equilibrium) statistical mechanics is to relate the microscopic prop- 
erties of a fluid (such as a knowledge of the interaction potential referred to above) to a set of 
coefficients which appear in the macroscopic equations of motion for the fluid, and describe the 
transport (motion in space) of conserved quantities; they are therefore known as transport coeffi- 
cients. These include coefficients describing diffusion (transport of mass), viscosity (transport of 
momentum) and heat conduction (transport of energy in the form of heat . 

1.1.3. Diffusion and Brownian motion 

In this thesis we discuss only diffusion in detail. This is one of the most fundamental phenomena 
in fluids and refers to the spreading out of matter which is initially confined to a small subregion 
of a system. This is the process by which concentration gradients in the system are smoothed out, 
and hence one of the ways in which a system returns to equilibrium. 

In seminal work published in his annus mirabilis 1905 IIEin85ll (see also IIGar85[ Chap. 1]), 
Einstein related diffusion to the motion of a large particle subjected to repeated impacts with 
molecules of the surrounding fluid. Regarding these impacts as random, due to the complexity 
of the presumed true microscopic dynamics of the huge number of fluid molecules, he showed 
that certain plausible assumptions led to an equation describing the probability distribution of the 
position of the particle; this equation is identical to the classical diffusion equation. The seemingly 
random motion of such particles had been observed in 1827 by Brown, a botanist, studying the 
motion of a pollen grain in water, and is hence called Brownian motion. For a recent experiment 
of relevance to this thesis see llGBF+981 . 

Einstein's derivation effectively models the motion of the particle as a stochastic process, 
namely a random walk; it is thus a mesoscale model intermediate between the true underlying 
microscopic dynamics and the observed macroscopic behaviour. 



These are the only transport coefficients for a simple iiuid, i.e. a iiuid consisting of one component. 
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1.2. Billiard models 

1.2.1. Toy models exhibiting diffusive properties 

Recently it has been realised that it is possible to study the statistical properties of some simple 
deterministic dynamical systems at the level of the full microscopic dynamics, and that these 
resemble to some extent those of diffusion. When this is the case, we talk about deterministic 
diffusion. As discussed extensively in Chap. |2l there is a hierarchy of related statistical properties 
which may 'look diffusive'. The concept of deterministic diffusion thus actually consists of several 
levels. 

Such systems can be regarded as 'toy' models to understand transport processes in more real- 
istic systems IIDor99ll . Examples include classes of uniformly hyperbolic one-dimensional (ID) 
maps (see e.g. IIKD99II and references therein) and multibaker models IIGas98l . Often, rigorous 
results are not available, but numerical results and analytical arguments indicate that diffusion 
occurs, for example in Hamiltonian systems such as the standard map IILL92I . 

1.2.2. Billiard models 

Billiard models constitute an important class of such models. Here, non-interacting point particles 
in free motion undergo elastic collisions with an array of fixed scatterers. Such models were 
introduced by Lorentz MLorOSL who modelled electron flow through an amorphous metal by point 
particles moving through a random array of hard spheres; such a model (with hard discs or spheres 
as the scatterers) is now termed a Lorentz gas. Several example trajectories of a periodic Lorentz 
gas are shown in Fig. 11.11 

Lorentz gases have been used to model neutron transport in dense media IICZ67II and can be 
viewed as modelling the flow of a dilute gas of light particles through a gas of heavy particles, 
in the limit where the ratio of masses of the light particles to those of the heavy particles goes to 
infinity IICC70II . The Lorentz gas is also a basic model in kinetic theory, since certain questions 
are simpler to answer in this context BvZvBDOOl IHau74l . 

1.2.3. Hard-sphere fluids as billiard models 

Furthermore, a hard-sphere fluid can be regarded as a bilUard model in a high-dimensional phase 
space. The simplest example of this is a 2-disc periodic fluid consisting of two discs on a torus, 
or equivalently two discs in each copy of a periodically-repeated unit cell, as shown on the left of 
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Figure 1.1: Sample trajectories in a periodic Lorentz gas. Each trajectory emanates from one point in the 
central unit cell, with only the initial velocities being different. The lattice of scatterers extends throughout 
space, but for clarity only a portion is shown. 
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Fig. 11.21 Call the discs A and B, with radii and rg. Suppose we stand at the centre of disc A 
and look at the motion of disc B bouncing off disc A. We see the same dynamics as we would if a 
point particle located at the centre of B were colliding with a fixed disc of radius + located at 
the centre of disc A. In this way the dynamics of the periodic fluid is equivalent to the dynamics 
of the periodic Lorentz gas model shown on the right of Fig. 11.21 



More generally, any hard-sphere fluid in a torus can be regarded as a billiard in a high- 
dimensional phase space, as follows. Consider N hard spheres in a periodic box QcM!^, with po- 
sitions q, E Q, velocities v,- G M'^, masses m; and radii r,. The vector (qi, . . . ,qA?,vi, . . . ,\n) G -M' 
then specifies the instantaneous state of the system, where M' := x is the phase space of 
the dynamical system. 

The total momentum is conserved, so we can change to a frame of reference in which the 
centre of mass is fixed at and the total momentum is 0, allowing us to restrict attention to a 
reduced phase space M of lower dimension ||Sza93l . Furthermore, regions of M corresponding 
to configurations where the hard spheres overlap are not allowed; these excluded regions are high- 
dimensional cylinders in A4. The dynamics of the hard-sphere fluid then corresponds exactly to 
free motion in together with elastic collisions on these cylinders, i.e. to billiard dynamics in 
the phase space M. IISza93ll . 
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1.2.4. Transport in billiards 

Billiards can be regarded as the simplest physical systems in which diffusion, understood as the 
large-scale transport of mass through the system, can occur: as pointed out in MBunOOII . all that 
is required for diffusion to be possible is some mass which can move through the system, which 
is exactly the situation in billiards. Other transport processes have also been studied in billiards, 
e.g. electrical conduction IICELS93b[|CELS93aL heat conduction IIAACG99II and viscosity IIBS961 
IBun00llVG03l . 

In this thesis we consider only diffusion in periodic billiard models, where the scatterers form 
a periodic array. In this case the dynamics in the extended (unfolded) system can be obtained by 
looking at the dynamics on a single unit cell with perodic boundaries, i.e. a torus with scatterers 
removed, and keeping track of which lattice cell particles are in. The region Q exterior to the 
scatterers on the torus is then called the billiard domain. Since the particles are non-interacting, 
it is usual to set all velocities to 1 by a geometrical rescaling, although in Sec. I4.7l we discuss the 
effect of a Gaussian velocity distribution. 

1.3. Classical diffusion equation 

Since we wish to describe dynamical systems as diffusive if their statistical behaviour looks like 
diffusion, we briefly review some features of the classical diffusion equation (also called the heat 
equation). This is a partial differential equation which models the empirically observed diffusion 
(flow / transport) of matter from regions of high concentration to regions of low concentration, and 
is one of the classical field equations of macroscopic physics. We begin by deriving the diffusion 
equation via a macroscopic balance equation. 

1.3.1. Balance equations 

Let p be the density (per unit volume) of an extensive quantity. Then p(r,?) : M'^ x ^ M is 
a scalar time-dependent field depending on position r G R'^ in space and time t, where d is the 
number of spatial dimensions. Examples of such fields are the mass density of a substance, the 
heat content, the concentration of a chemical species, and the x component of the momentum. We 
assume that p is sufficiently smooth that we can differentiate with respect to space and time. 

Consider a fixed region V (a 'volume') in the space W^, with outward unit normal vector n and 
boundary S :=dV. The total amount of substance (or in general of the extensive quantity of which 
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p is the density) inside V at time t is JyP (r, t) dr; this is the mass inside V if p is the mass density 
of a gas, for example. 

This amount of p inside V can change over time by exactly two (local) mechanisms: p can 
be created or destroyed inside V, with source strength a, the amount created per unit volume, per 
unit time; or p can flow across the boundary of the volume V, with flux vector J per unit area, per 
unit time. This neglects non-local effects such as radiation. 



Thus 

d 
df 



I p(r,Odr= I a{r,t)dr- I JndS. (1.1) 

V V dv 

Since V is fixed, we can move the time derivative inside the integral. The divergence theorem then 
gives 



I ^(r,Odr= I a(r,Odr- J V-Jdr, (1.2) 



dp 
dt 

V V 





for any fixed volume V. This implies thai 

|^(r,0 = a(r,f)-V-J(r,f). (1-3) 

The equation (11.31 ) is termed a balance equation IIKP98II . It forms the basis for deriving macro- 
scopic evolution equations for any scalar field, and hence also vector and tensor fields by consid- 
ering components. In order to derive such equations, we must specify a and J in terms of other 
known quantities using constitutive equations, which model the 'constitution' (behaviour) of a 
substance, to get a closed system of equations: see e.g. IIKP98II . 

1.3.2. Derivation of the (anisotropic) diffusion equation 

We now specialise to the case where p is the concentration, or mass density, of a substance which 
is diffusing. If the mass of the substance is (locally) conserved, we have a = oE We must now 
specify J, a question which is considered at length in non-equilibrium thermodynamics lldGM84l . 



"^Taking all terms to one side gives an equation of the form Jy /(r,r)dr = 0. If / is not (almost) everywhere 0, 
then it must be positive (without loss of generality) on some region V+ of non-zero volume. Then /y+ f(r,t) dr > 0, 
contradicting the assumption that the integral equal zero over any volume. 

^This is not the case, for example, if we have several reacting chemical species. The source Oa of species A then 
contains terms describing the creation and destruction of A in the reactions (the total mass, however, being conserved). 
We then obtain a set of reaction-diffusion equations describing the spatio-temporal evolution of the concentration fields. 
Such equations can yield interesting spatial patterns [KP98J. 
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We consider only the (simplest) case, when Pick's law is obeyed, so that the flux J is a linear 
function of the concentration gradient, Vp . (This can be viewed as a first-order approximation un- 
der the assumption of local thermodynamic equilibrium and slowly varying concentration profile 
lldGM84llKP98]| .) Hence 

J = -D-Vp, (1.4) 

where D is a second rank tensor called the diffusion tensor. The minus sign accounts for the em- 
pirical fact that matter diffuses from regions of high concentration to regions of low concentration. 
This is an example of a constitutive equation. 

Assuming that D is independent of position, we have 

V . J = -Y^bijdijp = -Y^DijdijP, (1.5) 

i,j hi 

where D,y are the components of the tensor D with respect to a Cartesian coordinate system, and 
dip := We have defined the (symmetric) diffusion tensor D as the symmetric part of D, i.e. 
D = ^(D + D^), with components Dij. The antisymmetric part does not play any role, since the 
matrix of second partial derivatives of p is symmetric. Substituting in (11.31) gives the (anisotropic) 
diffusion equation 

^=LDijd^jp. (1.6) 
We can also write this independently of coordinate system as 

|^ = V.(D.Vp). (1.7) 



1.3.3. Solution of the diffusion equation in an unbounded domain 



We consider the diffusion equation in one dimension, with D independent of space: 



We solve equation dl.SI ). a linear partial differential equation with constant coefficients, using 
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a standard Green function method. For G R we define tlie Fourier transform of p at time t by 

oo 

p(t;k):= J e-'^'p{t;x)dx. (1.9) 

Fourier transforming dl.SI ) gives 



5? 



(1.10) 



an ordinary differential equation for p{t;k), the solution of which is 



-Dk't 



p{t;k)= p{0;k)e 
Taking the inverse transform of this product gives a convolution: 

oo 

p,{x) = p{t-x) = ^ I e+'^'p{t;k)dk 

k— — oo 

= [po*G']{x), 



(1.11) 



(1.12) 
(1.13) 



where po(x) := p(0;x) is the initial concentration distribution and G' is the Green function (or 
propagator) for the diffusion equation on an unbounded 1 -dimensional domain. This Green func- 
tion is given by the Gaussian 



G'{x) ■.= G{t;x) :-- 



1 



yjAnDt 



exp 



ADt 



with mean and variance 



Var[G'] := j x^ G' {x)Ax = lDt 



at time t, and Fourier transform 



G{t;k)=& 



-Dkh 



(1.14) 



(1.15) 



(1.16) 



We recall that the convolution operation is defined by 

oo 

{u*v){x):= / u{x-y)v{y)dy. 



(1.17) 
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Multidimensional diffusion equation In the multidimensional case, since D is a symmetric 
tensor it is possible to choose an orthogonal coordinate system in which D is represented by 
a diagonal matrix. In these new coordinates, the solution is a product of solutions of the ID 
equation. Reverting to the original coordinates gives a multi-dimensional Gaussian IIGS92II for a 
Dirac-delta initial condition. 

1.3.4. Solutions of the diffusion equation as probability densities 

Henceforth we regard the diffusion equation as describing the evolution of probability density 
functions representing probability distributions, as follows. 

Let the initial condition at time ? = be po. Physically we are interested in non-negative po 
with finite mass, i.e. 

oo 

J po{y)dy<^. (1.18) 

y— — oo 

By normalising if necessary we can instead assume that 



I po{y)dy = l. (1.19) 



It then follows that J'^„Pt{y)dy = 1 and that the solution remains non-negative for all times t, so 
that we can regard the diffusion equation as describing the time evolution of probability densities. 

1.3.5. Calculation of moments from Fourier transform 

We can calculate moments of a distribution directly from its Fourier transform as follows (see e.g. 
llRD77llBal97l ). 



Differentiate ( 11.91 ) with respect to k to get 

o 

dp{t;k) 



dk 

and hence 

dp{t;k) 



k=Q 



j -ix&-'^''p{t;x)dx, (1.20) 



J -ixp{t;x)dx=: -i{x)t. (1-21) 



dk 

Here we denote by {f{x))t the mean of the function f{x) with respect to the distribution at time t. 
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Similarly, higher moments can be obtained by differentiating repeatedly with respect to k: 

oo 

d"'p{t;k) 



dk" 



-i)™x™ p{t;x)dx= (-i)™ {x^)t. (1.22) 



k=Q 

— OO 



In a similar way, in the full multi-dimensional, anisotropic diffusion equation (11.61) for p{t;x) 
in n dimensions, we use the multi-dimensional Fourier transform 



oo CO 



p{t;k):= I ■■■ J e-'^''p{t;x)dxi---dxn, (1.23) 



— oo — oo 



where k • x := ^"^j kjXi. Then we have, for example, 

52p(f;k) 



{XiXj)t 



(1.24) 



k=0 



dkj dkj 

with higher moments calculated in an analogous way. 

1.3.6. Asymptotic behaviour of moments of solutions of the diffusion equation 

We assume that po decays sufficiently fast at infinity for the relevant integrals to exist. Physically 
relevant sufficient conditions for this are, for example, that the initial condition has compact sup- 
port, i.e. it vanishes outside a finite interval, or that the initial condition is exponentially localised, 
in the sense that 

Po(j)<e-^W (1.25) 

for some constant K >0. 

In this section we denote partial differentiation with respect to ^ by a subscript k, so that 
pk :=dp/dk. 



First moment We first calculate the time dependence of the first moment (i.e. the centre of 
mass) of the solution of the diffusion equation. Differentiating (11.111) (the solution of the diffusion 
equation expressed in Fourier transforms) once with respect to k, we have 



Pk{t;k) = {-2kDtpo{k)+pok{k)} , 



(1.26) 
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SO that using equation (11.21b twice, we obtain 

{x), = ipkit;0) = ipo^(O) = ip;t(0;0) = (^)o. (1.27) 

Hence the mean position (centre of mass) is constant. So the second moment is the first non-trivial 
one. 

Second moment To simplify the calculation of the time evolution of the second moment, we 
now change to a coordinate system with the origin located at the centre of mass, i.e. we assume 
without loss of generality that (x)o = 0. 

Differentiating (11.111 ) again with respect to k gives the second derivative 

Mt-,k)=c-'''''{po,,{k) + 2po,{k).{-2kDt)+po{k) [-2Dt + {-2kDt)^]}, (1.28) 
so that 

{x^), = -Pkk{t\0) = -po^^.(0) +2Dfpo(0) = {x\ + 2Dt. (1.29) 

(Note that po(0) = / 1^ '^^ Pq{x)Ax =1.) Since po was assumed to have compact support, po is 
infinitely differentiable, by regularity results for Fourier transforms IIKat04l . Hence {x^)t is a 
straight line with slope 2D, which does not pass through the origin unless the initial condition is a 
Dirac delta (since the initial variance is only if the initial distribution is concentrated on a single 
point); the variance hence grows asymptotically linearly. 

From the above result we can derive a relation between the diffusion coefficient and the rate of 
growth of the variance. Dividing ( 11.291 ) by t and taking the limit gives 

D = lim^. (1.30) 
It 

This result is known as the Einstein relation, since it was first obtained by Einstein IIEin85ll . We 
also have 

D = \im~y)t, (1.31) 
r^cxj 2 at 

but note that the first limit can exist when the second one does not, for example if {x^)t = sin t. 
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Similarly in the multi-dimensional case we find 



If the system has enough symmetry (Sec. 13.61 ). then we have D„- = D for all /, with all other 
components equal to 0, and 

D = \mi\£, (1.33) 
2dt 

where d is the spatial dimension and := x • x = I^f^iX?. We remark that this reduction of the 
diffusion tensor to a multiple of the identity tensor also occurs if the system is isotropic, i.e. has 
the same properties in any direction in space, but isotropy is a stronger condition than necessary 
for this reduction to happen. 



1.3.7. Convergence of solutions of the diffusion equation to Gaussians 

If the initial condition for the diffusion equation is po at time ? = 0, then the solution at time t is 
given by the convolution 

oo 

p,(x)= j po(3^)-^e-(^->')'/4^M);. (1.34) 

y— — oo 

We are interested in the shape of the distribution for long times. In App. |A] we show that 
for suitable initial data the solution converges to a Gaussian when appropriately rescaled. This 
convergence requires the limiting function to be non-degenerate. Since pt tends to pointwise 
as f ^ oo, we must first rescale p,. We know from the above argument or from dimensional 
considerations that ~ T , where X is a typical lengthscale and T is a timescale. Hence we 
rescale x by ^/t, putting 

p,{x) := p{t\x) := ^ft .p{t\x^ft). (1.35) 
The first factor of ^/t is to normalise the integral of Pf to 1 . 

In Chap. |4] we study the convergence to a limiting normal distribution of rescaled distribution 
functions in the context of billiards. For comparison, in App. |A]we consider the convergence of 
rescaled solutions of the diffusion equation to the hmiting Gaussian. We show that Pt {x) converges 
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pointwise to a Gaussian with variance ID: 

pt{x)'^ g2D{x), (1.36) 

where 

is the Gaussian density with mean and variance 2D. Furthermore the convergence is uniform 
with 

\Pt{x)-g2D{x)\^j, (1.38) 
for all X and some constant C which is independent of x. 

1.3.8. Numerical confirmation of 0{t^^ ) decay 

We numerically confirm the 0{t^^) decay found in App.lAlby considering particular initial distri- 
butions. The first initial condition we consider is 

Pi{x):=\8{x-2) + l8{x+\), (1.39) 

where 5{x — xa) is the Dirac delta function at position xq. The diffusion equation with this initial 
condition is analytically soluble, with solution 

Pt{x) = \Gt{x-2) + lGt{x+\), (1.40) 

so that the solution is smooth and rapidly decaying for any t > and hence fits into the class for 
which we can prove the 0(f^^) convergence. 

To verify the convergence numerically, we calculate ||p, — §20! I = sup_f£]g \pt{x) —gioix)], set- 
ting Z) = 1 by rescaling time. Plotting the logarithm of this distance against the logarithm of time 
gives a straight hne with slope —1, confirming the t^^ decay, as shown in Fig. 11.31 

Another initial condition for which the equation is exactly soluble is 

P2ix) = ^i[-K,K]i^), (1.41) 
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Figure 1.3: Uniform distance £, (defined in App.|A]i from the limiting Gaussian of the rescaled solution of 
the diffusion equation with initial conditions pi(x) = ^G,(x — 2) + |G,(x+ 1) and P2{x) = j^l[_K^K]ix)- 



with solution IICra75ll 



2y/DtJ \2y/DtJ 



(1.42) 



where 



X 

erf(x):=^ U-'^'l^dX 



(1.43) 



is a special function called the error function. Again the numerical 0{t decay rate, as shown 
in Fig. 11.31 confirms the analytical result. 

1.3.9. Hallmarks of diffusion 

We conclude from the above remarks that hallmarks of diffusion are the following features of the 
asymptotic behaviour as f ^ oo; 

• mean squared displacement growing asymptotically linearly in time t, with constant of pro- 
portionality 2D, where D is the diffusion coefficient; and 



convergence of the -i/F-rescaled position distribution to a non-degenerate Gaussian. 
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Here by 'non-degenerate' we mean that the Umiting Gaussian has non-zero variance, or in higher 
dimensions that the covariance matrix is positive definite. 

These are the statistical features which we will use to characterise a process as diffusive, to- 
gether with a third dealing with convergence of rescaled paths to Brownian motion. 

We remark that in general inhomogeneous physical systems, the rate of diffusion characterised 
by the diffusion coefficient can vary over the system. In such a case, we can say that we have 
diffusion if Pick's law holds locally, i.e. if the flux of the diffusing quantity is locally proportional 
to its gradient IICra75l . 

1.4. Statistical properties 

Many dynamical systems, including important classes of billiards, are chaotic, in the sense that 
trajectories emanating from nearby initial conditions separate fast as time evolves. Any physical 
measurement has a limited precision, so that if we repeat an experiment several times then we 
cannot be sure that the initial condition is identical each time, but only that it lies within a certain 
tolerance. Averaging over the results of the experiments thus corresponds to averaging over a 
small ball of initial conditions as it spreads out over possible states of the system. This argument 
motivates the need for introducing probabilistic notions, namely ergodic and statistical properties, 
to study deterministic dynamical systems which are chaotic; see Chap.|2] 

One of the main reasons for the interest in billiard models is the possibility of obtaining rigor- 
ous results on their ergodic and statistical properties: see e.g. IIBS811|BSC91|| . The techniques are 
most highly developed for dispersing billiards such as the Lorentz gas, where the curved, convex 
scatterers cause nearby trajectories to separate exponentially fast. In fact, Lyapunov exponents 
(which measure the rate of separation) exist and are non-zero almost everywhere for the billiard 
map, so that the system is hyperbolic (chaotic). Due to the Hamiltonian nature of the system, 
the Lyapunov exponents come in positive-negative pairs, so that at least one Lyapunov exponent 
is positive almost everywhere IICMOli [Gas98ll ; further, the Kolmogorov-Sinai entropy, which 
measures the rate of generation of information in time, is positive IIGas98ll . These are standard 
indicators of the chaotic nature of the system IIER85II . 

Hard-sphere fluids are only semi-dispersive, due to the flat, neutral directions along the cylin- 
der axis in the many-particle phase space; this makes their rigorous analysis much harder: see 
e.g. IIBal99ll for a review. Nonetheless, a rigorous proof of the celebrated Boltzmann ergodic hy- 



1.5. Overview of thesis 



17 



pothesis has recently been achieved for hard discs {d = 2) and spheres {d = 3): see references in 
Sec. [2311 

The main focus of this thesis is deterministic diffusion, a statistical property of certain dy- 
namical systems, including certain classes of billiards. A definition often used in the physical 
literature is that a system is diffusive if the mean squared displacement grows proportionally to 
time t, asymptotically as f — > oo. However, there are stronger properties which are also charac- 
teristic of diffusion, which a given system may or may not possess: (i) a central limit theorem 
may be satisfied, i.e. rescaled distributions converge to a normal distribution as f ^ oo; and (ii) the 
rescaled dynamics may 'look like' Brownian motion. See Chap.[2]for details. 

Recently there has been much interest in the question of which microscopic features are neces- 
sary for a system to exhibit strong ergodic and statistical properties. The proofs of these properties 
for dispersive billiards depend crucially on the fact that they are hyperbolic, but numerical evi- 
dence has been given that systems with weaker chaotic properties may also show strong statistical 
properties, for example the polygonal billiard channels studied in Chap. [5] 

1.5. Overview of thesis 

Chapter|2]reviews concepts from dynamical systems, probability theory and ergodic theory which 
we require to discuss deterministic diffusion. Chapters [3]-[6] present our results on deterministic 
diffusion in three types of billiard model: a 2D periodic Lorentz gas with a two-dimensional 
parameter space (Chaps. [3]and|4l); several classes of polygonal billiard channel (Chap. [5ll; and 
3D periodic Lorentz gases (Chap. H]). Conclusions are presented in Chap. |71 followed by several 
appendices giving technical details of results needed in the main text; we expect that many of these 
results are already known, but in several cases we were unable to find suitable references. 

The discussion in each chapter can be thought of in terms of three related threads. The first is 
mathematical: to what extent do rigorous results hold beyond their immediate range of applicabil- 
ity, and can we understand more precisely the statistical behaviour described by those results? The 
second is physical: how do geometrical features affect the dynamics and statistical properties of 
the system? The third is statistical: how can we best estimate statistical properties from numerical 
data? 



Chapter 2 



Statistical properties of dynamical systems 



In this chapter we review in what sense deterministic dynamical systems 'look Uke' stochastic 
processes. We make a (somewhat arbitrary) distinction between ergodic properties and statistical 
properties: the former deal with general theorems referring to large classes of observables, whereas 
the latter hold only for observables with a certain degree of smoothness, are related to rates of 
convergence in ergodic theorems, and are of physical relevance. 

2.1. Dynamical systems and stochastic processes 
2.1.1. Invariant measures 

Let ^ : M Mhe the flow of a continuous (t G M) or discrete (t G Z) dynamical system with 
phase space M. Measures fi are defined on some (7-algebra B on M. Our dynamical systems 
will be defined on metric spaces, so that we can take, for example, B to be the Borel measurable 
subsets of M. Then any physically relevant subset of M will be in B, and from now on we will 
usually refer to subsets of M without expUcitly mentioning B. 

Suppose that we start the evolution of the dynamical system with a distribution of initial con- 
ditions described by the measure Hq. This distribution will evolve in time to the measure [Xt given 
by 

^t:=m4^y, M^(A):=/io(<J>"'(A)); (2.1) 

here the right hand side defines the meaning of the push-forward {^)*. 

A particularly simple and interesting case occurs if the measure is preserved by the system, 
or is invariant, as follows. We say that jx : = Hqis invariant with respect to the flow if, for all t, 
(<I>')*(/i) = [X, or equivalently if 

H i^-' (A) ) = At (A) , for all A G (2.2) 

Here we also require that ^ be measurable with respect to B, i.e. such that 0~'(A) G B for all 
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A G B. In the discrete time case, where the dynamics is given by a map T : M ^ M, it is enough 
to have |i (r^^ (A)) = /^(A) for all measurable sets A. 

Natural invariant measures Some systems possess natural invariant measures if the dynamics 
preserves some structure. The main class of interest to us is Hamiltonian systems, which preserve 
Liouville measure: see Sec. l2.5l in the case of billiards and IIKH95[ Chap. 5] for other examples. 
We are then most interested in statistical properties with respect to these measures. 

2.1.2. Dynamical systems with invariant measures as stationary stochastic processes 

Fix an observable f : M. ^M. (that is, a quantity that we could in principle measure when the 
system is in different states in phase space) and look at X, := /o : ^ M. Given a measure 
pL := /^o on Al, we can regard (Xf ) as a collection of random variables indexed by time, so that we 
have a stochastic process IIGS921I . 

If the measure /i is invariant, then the process (Xf) is stationary, which means that for all n, all 
ti,...,tn and all h>0, the families 

(X,,,...,X,„) and (X,,+,„...,X,,+,,) (2.3) 

have the same joint distribution as each other IIGS92II . The proof is as follows. For two times t and 
s and two sets A,B G ;B, we have 

P {fo^'eAJo^^ G B) = Ai {^-'{r\A))n^-\f-\B))) , (2.4) 
which by the measure-preserving property of <I> is equal to 

= M (<!>-'■ {^-'{f-\A))n^-'{f-\B)))) = P {fo^'+'- e A,fo^^+'- G B) . (2.5) 

An analogous argument holds for n times and n sets; hence the random variables {X, = f o^*) 
form a stationary stochastic process. 

2.2. Ergodic properties 

Ergodic properties with respect to an invariant measure can be thought of loosely as a measure- 
theoretic description of chaoticity. There is a hierarchy of increasingly strong properties; the ideas 
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originated in statistical mechanics: see e.g. ||Dor99ll . A general reference is, for example, IICFS82I . 

2.2.1. Ergodicity 

In the billiard systems we study there is a natural invariant measure /i, namely Liouville measure 
(see Sec. 12.51 ). which we regard as fixed in the following. 

Suppose that the invariant measure jJ. is finite, i.e. jU(A^) < °°. Then we can normalise to get 
a probability measure with /x(A^) = 1. The flow <!>' is ergodic with respect to ji ij^for all A C A^, 
we have 

^~'{A)=A pl{A)=0 ox pl{A) = \, (2.6) 

i.e. any set A which is invariant under the flow has measure 0, so that it is trivial from the point 
of view of measure theory, or has measure 1 , so that it covers the whole space except for a set of 
measure 0. 

2.2.2. Birkhoff ergodic theorem 

Let fJL be an invariant probability measure for the dynamical system O', and let / : ^ M be 
integrable. Then Birkhoff 's ergodic theorem states (see e.g. IICFS82II ) that the limit 

T 

/(^):= hmi f{fo^%x)dt (2.7) 



exists for almost aWx^ M with respect to the measure }X, and we call /(x) the time average of /. 
Furthermore, if also the system is ergodic, then 

f{x) = if)^ := J /dM a.e., (2.8) 
M 

SO that time averages are almost everywhere (a.e.) equal to the space average (/) with respect to 
the ergodic invariant measure /i. 

The motivation for this theorem came originally from Boltzmann, who had the idea that a 
sufficiently complicated dynamical system should explore the whole accessible phase space, an 
idea known as the ergodic hypothesis. This has recently been proved for hard-sphere fluids in a 
series of papers ISiml ISim04[ ISim031 [SS99ll . 

' We continue to suppress the necessity of iiaving A e S. 
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2.2.3. Mixing 

The flow ^ is mixing with respect to the invariant measure jJ. if for any sets A,B C A^, we have 

^iAn^'iB))'^^iA)^iB). (2.9) 

Mixing implies ergodicity. An important interpretation for us is that mixing is equivalent to weak 
convergence of densities: see App.lCl 

There is also a notion of weak-mixing, intermediate in strength between ergodicity and mixing. 
For discrete time systems we say that a transformation T preserving a measure /i is weak-mixing 
if for any A,B C we have 

lim - y \pi{T-'Ar\B) - ii{A)ii{B)\ =0. (2.10) 

For comparison, there is a characterisation of ergodicity in similar terms, stating that T is ergodic 
if for any A,B C M we have 

lim-yn{T-'Ar\B)=n{A)n{B). (2.11) 

See IIWal821 ICFS82II for detailed comparisons of the different notions, and IICFS82II for formula- 
tions for continuous-time systems. 

2.2.4. K-systems 

K-systems have very strong ergodic properties: in particular their Kolmogorov-Sinai entropy is 
positive. See e.g. llSin94[ p. 71] for the definition. K-systems are multiply mixing (a generalisation 
of mixing), mixing and ergodic. 

2.3. Statistical properties: probabilistic limit theorems 

We now turn to statistical properties, reviewing the application of results from the theory of sta- 
tionary stochastic processes to the context of dynamical systems. 

Diffusion in billiards concerns the statistical behaviour of the particle positions. Denoting the 
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position at time t by x, and restricting attention to the first component X; , we can write 

t t 
xt = jvi{s)As + XQ = j fo^\.)ds + XQ, (2.12) 



where / = vi, the first velocity component. This expresses x, solely in terms of functions defined 
on the torus, so that in a sense we have reduced a spatially extended problem (spreading out 
over an infinite lattice) to a problem on the torus. Equation (12.121 ) shows that the displacement 
Axf :=Xt —xqis in some sense a more natural observable than the position Xt in this context. 

In the above we are regarding vi : — > R as the observable 

vi(«)=vi(^,v)=vi (2.13) 

which returns the first component of the velocity of the initial condition, so that vi o^[(o) is the 
first velocity component at time t starting from the initial condition CO = {q,v) € Ai. 

More generally, we can consider integrals of the form in ( 12.121 ) over other observables / : 
7W ^ M; these are important in the study of other transport processes, for example HBunOOi The 
question we wish to answer then concerns the distribution of accumulation functions of the form 
lICYOOt 

t 

Sr{-):= lfo^\-)ds, (2.14) 


in particular in the limit as ? ^ oo. 

The integral in (12.141 ) is a continuous-time version of a Birkhoff sum 

Sn=f + foT + foT^ + ...+foT"-' =Xi + ---+Xn, (2.15) 

where the X„ := /o T"^^ are stationary random variables. We are thus interested in statistical prop- 
erties, for example means and shapes of limiting distributions, of sums and integrals of stationary 
random variables. 

Intuitively, if the correlations between the Z, decay sufficiently fast, then they are asymptoti- 
cally independent, and we can hope that the classical hmit theorems for independent and identi- 
cally distributed random variables can be extended to the stationary case. The Bernstein method, 
summarised below, is a rigorous justification for this. 
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2.3.1. Averages 

The simplest statistical propeities to study are averages of observables. If / : ^ M is an 
observable on the phase space of the dynamical system, then we denote the mean of the observable 
/by 

(/) :=E^[/] = |/dAt = 1 f{co)d^{co). (2.16) 

M M 

Here, jJ. is the distribution of initial conditions co € Ai and [/] denotes the expectation of the 
random variable / with respect to the probability measure jJ.. 

We will be interested in the evolution over time of such averages when the observable / in- 
volves the flow <I>^ A key role is played by the mean squared displacement at time t, denoted 




{Ax\:={{Axtf)^ = {{x,-xof)^ = (lJv,o^^{.)ds\ \ . (2.17) 



In the final expression in (I2.17I ). the dot denotes the variable O) G over which the average (•)^ 
is taken. We will usually think of this mean squared displacement as a function of time t, so that 
it is convenient to use the notation (Ac^),, which makes this dependence explicit. 

The physical interpretation of the above definition is as an average over initial conditions of the 
time-dependent observable. An alternative point of view is to regard {Ax^)t as an average over the 
evolved probability distribution /i, at time t starting from a distribution /Iq at time 0; the average 
is then over a fixed observable at points determined by evolving the initial conditions in time. 

2.3.2. Central limit theorem: independent, identically distributed case 

Let (Xi) be independent and identically distributed (i.i.d.) random variables with mean and 
variance < o°. We are interested in statistical properties of the accumulation function 

S„:=Xi + ---+X„, (2.18) 

as n ^ oo. We have 

E [Sn] = 0; Var [S,,] = « Var [X] , (2.19) 
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SO that if we normalise by yjn, setting 5„ := Sn/^/n, then we get 

Var[5J 



Var [Sn] = Var 



a^. (2.20) 



Since the variances of the iS^ are now independent of ?i, we may hope that the distributions converge 
in shape to some limiting distribution. 

With the above conditions, the classical central limit theorem (see e.g. IIFel711 IDur96ll ) states 
that 

— oo 

so that the rescaled accumulation functions converge in distribution to a normal distribution. 

We recall the definition of this notion of convergence IIGS921lBil68ll . Let X„ be a sequence of 
real- valued random variables with distribution functions F„, so that F„{x) := P {X„ ^ x). We say 



that the sequence X„ converges in distribution to the random variable X, with distribution function 
F{x) := P (X ^ x), written X„ — > X, if F„{x) —>■ F{x) at all x where F{x) is continuous. There is 
a generahsation of this concept which applies in much more general situations: see App.lDl 



2.3.3. Central limit theorem: stationary case 

If now the X,, are no longer independent but they are stationary, we have 

n-l 

Var[5„] = «C(0) + 2 £ in-j)CU) (2.22) 

where the autocorrelation function of / is 

C{n) :=E[XoX„] = {f.{foT"))-{f)\ (2.23) 
and (•) = E [•] both denote averages (expectations) over the invariant measure. It follows that 

^^"-^a':=m+2tcU), (2.24) 



provided jCj exists. A sufficient condition for thij^ is that Cy = o (^j ^ for some £ > 0. 



Note that Cj = o (7 ^) is not sufficient, since Y,j jCj is divergent for Cj = j ^. On the other hand, Cj = o iyj ^ ^) 



is not necessary either, since e.g. for Cj = j ^(logy) ^ the sum Y-jjCj is convergent for p > 1 IRud76l . 
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Figure 2.1: Partition for the Bernstein method. 



Now we may again hope/expect that 




(2.25) 



We can also regard the central limit theorem as describing the distribution of 'fluctuations' around 
the mean in the Birkhoff ergodic theorem. (Note that the Birkhoff ergodic theorem is a version of 
the strong law of large numbers for stationary processes IIDur96ll .) 

2.3.4. Central limit theorem for stationary processes: Bernstein method 

The Bernstein method can be used to prove results of central Umit theorem (CLT) type for station- 
ary processes, based on the idea that if E [XqX,,] decays sufficiently fast then the X, are asymptoti- 
cally independent and we can reduce to the case of independent random variables. 

The method works as follows (see e.g. ||Dur96[ Chap. 7] and IIB 116811 '). Split the sum 5„ up into 

alternate blocks t,j of length p and r]j of length q, so that t,\ =X\-\ \-Xp, r]\ = Xp+\ H f- 

t,2 = Xp+^+i + • ■■X2p+q, etc.; see Fig.O 

If q is sufficiently large then the i§y are almost independent; on the other hand, if q is small 
compared to p then the sum Y,j Vj of the small blocks is small compared to the total size of the 
sum 5„. In a suitable limit we can thus reduce to the independent case, with an error which tends 
to 0. 

2.3.5. Rate of mixing 

For the above to work, we need to show for a given stationary sequence that the correlations C„ 
decay fast enough, i.e. that we have a fast enough rate of mixing. There are various sufficient 
'mixing conditions' for the central limit theorem to hold: see e.g. IIDur961 lBil681 IDav94ll in the 
context of probabihty theory, and IIDen89[|Liv96ll in the specific context of dynamical systems. 
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2.3.6. Functional central limit theorem 

The central limit theorem quoted above deals with one-dimensional distributions. There is a multi- 
dimensional version, where multi-dimensional distributions converge to multi-dimensional normal 
distributions. If also a further technical condition, called tightness, is satisfied, then the probability 
distribution on the space of continuous paths which is induced by the dynamical system converges 
weakly to a continuous-time stochastic process. If further the correlations behave correctly, then 
the limit will be the particular case of Brownian motion. Such a result is known as afunctional 
central limit theorem (i.e. a version of the central limit theorem for paths represented by functions). 
This is also called the weak invariance principle. 

For further details in the context of diffusion see Sec. 12.41 and App.lDl 

2.3.7. Almost-sure invariance principle 

The strongest type of statistical property is known as an almost-sure invariance principle, since 
it says that a stochastic process can be written, almost surely, as the sum of a Brownian motion 
on a suitable space, together with an error which can be bounded in a precise way: see e.g. 
llPS75llDP84llMT04ll . 

2.3.8. Transition from discrete-time to continuous-time 

In billiards, we are interested in statistical properties of the physical continuous-time dynamics. 
Currently these are not directly available, due to the lack of information on decay of correlations 
in continuous time; rather they are derived from the results on the billiard map, using the fact that 
the flow is a suspension flow over that map. 

We refer to App.lBjfor the definition of suspension flows and for a statement of the key theorem 
relating the central limit theorem for the flow to that of the map. We give here a derivation of the 
Green-Kubo formula from the Einstein formula for the diffusion coefficient in continuous time, 
valid only if the velocity autocorrelation function decays sufficiently fast, which for billiards has 
not yet been proved BCYOOI . References containing similar derivations include llvB82[ iDetOOl 
IBY911IRD77I . We remark that recent result have shown the exponential decay of correlations 
in continuous time for Holder observables of Anosov geodesic flows l|Dol98[ IChe98i and contact 
flows BLivL 

^ [Stretched exponential bounds on correlations for Holder-continuous observables for the 2D periodic Lorentz gas 
were recently proved in |Che07| .1 
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We take as a starting point the Einstein formula 

D=liml(Ax2),. (2.26) 

t^°° It 

Using Ax(?) =x{t) — x(?o) = ft'=Qv{t')At' , we have 

t t t t' 

{Ax^)t= J J {v{t')v{t"))dt'dt" = J J (v(T)v(0))d?'dT. (2.27) 
f'=Of"=0 f'=OT=f'-f 

Here we have used the fact that the averages are invariant with respect to time translation, due 
to the stationarity of the stochastic process; we also performed a change of variables from t" to 
z:=t'-t". 

We now change the order of integration, obtaining 

f+T t t 

{Ax^)t= J J (v(T)v(0))df'dT+ J J (v(T)v(0))df'dT. (2.28) 

t=-tt'=0 T=Of'=T 

Hence 

t 

{Ax^)t = 2 j {t-T:){v{z)v{0))dT:. (2.29) 

We can now define a finite-time diffusion coefficient D{t), an estimate for the (infinite-time) 
diffusion coefficient based on the information available up to time t. There are two possible ways 
to do this. Using the Einstein definition, we have 



(2.30) 



where C(t) := (v(t) v(0)) is the velocity autocorrelation function. Hence the diffusion coefficient 



D := ]imD(t) = [ C(i:)dT (2.31) 

f-»oo J 



exists if: (i) J^^QC{T)dz < °o, i.e. C is integrable; and (ii) j Jl^QzC{z)dT — >^ as ? — >^ oo. In fact 
(ii) follows from (i) by integrating by parts, so that a necessary and sufficient condition for the 
existence of the diffusion coefficient is that C{t) is integrable. 
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The equation (12.311 ). relating the time integral of the velocity autocorrelation function to the 
transport coefficient D, is a Green-Kubo relation ||Dor99ll ; Green-Kubo formulae can be found 
for all transport coefficients IIEM901 lGas98l . e.g. via linear response theory IIEM90II . expressing 
the transport coefficient as a time integral of the autocorrelation function of the flux of the quantity 
being transported. 



Another possible definition of a finite-time diffusion coefficient D, which is more closely re- 
lated to the numerical method we shall use later ||Bal971lAT87L is to define D using the local slope 
of the mean squared displacement (Ax^)f: 



From ( 12.291 ). we have 



{Ax^), = 2t J C(T)dT-2 J TC(T)dT, 



T=0 T=0 

so that the fundamental theorem of calculus gives 



(2.32) 



(2.33) 



D{t)= J C{T)dT + tC{t)-tC{t)= J C{r)dT. (2.34) 

T=0 T=0 



This definition avoids the slowly-decaying l/t tail of the first definition ||Bal97i . Further, if we 
have a bound on the rate of decay of correlations then we can use (12.341 ) to estimate the error 
|D — D(?)| resulting from using the finite-time diffusion coefficient D{t) instead of its infinite- 
time limit D. For example, if correlations decay exponentially as 



\C(t)\<Ke- 



(2.35) 



then we can estimate 



\D-D{t)\ 



C(T)dT 



■Z=t 



K 
— ( 
a 



(2.36) 



Hence this difference tends very rapidly to zero, as seen in numerical simulations. However, 
obtaining the constants K and a, even numerically, is difficult or impossible. 
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If, on the other hand, correlations decay only algebraically as 



-l-e 



(2.37) 



then 



\D-D{t)\ ^ 



C{s)ds 



(2.38) 



This error tends to zero very slowly for small e, so that the finite-time estimation method for any 
practical time has an additional uncertainty built in. 



2.4. Definition of deterministic diffusion 

There is no single definition of deterministic diffusion; instead, we say that a dynamical system is 
diffusive if one or more of the following hierarchy of statistical properties holds. The idea is that 
the statistical properties behave to some extent like those of solutions of the diffusion equation 
described in Sec. 11.31 or those of Brownian motion (see below). We use terminology suitable 
for diffusion in billiards, although these properties are also relevant for more abstract dynamical 
systems. 



(a) the mean squared displacement grows Unearly asymptotically; 

(b) the position or displacement distribution, rescaled by ^/t, converges in distribution to a 
non-degenerate Gaussian, i.e. a central limit theorem holds; and 

(c) the whole process, suitably rescaled, converges in distribution to a Wiener process. 



Properties (a) and (b) are based on the corresponding properties of the diffusion equation, 
although property (b) will in general only hold at the level of weak convergence, whereas in the 
classical diffusion equation there is pointwise (and even uniform) convergence of the densities. 
Property (c) says that the process, when suitably rescaled, looks like Brownian motion, which can 
be thought of as a probabilistic model of diffusion. We have the implication^ (c) =^ (b) =^ (a), so 
that (c) is the strongest property. We now give precise versions of these statements. 



'*[In fact, as was pointed out to me by Ian Melbourne, (b) does not imply (a) in general. However, in billiards, they 
tend to go together - to know which scaling factor to use in the central limit theorem, it is necessary to calculate the 
mean squared displacement.] 
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(a) Asymptotic linearity of mean squared displacement The limit 

2D:=lim-{Ax^)t (2.39) 

exists, so that the mean squared displacement (Ax^), := ([Ax(f)]^) (the variance of the displace- 
ment distribution) grows asymptotically linearly in time: 

{^x^)t ~ 2Dt as f ^ oo, (2.40) 

where D is the diffusion coefficient. Ind ^ 2 dimensions, letting Axi{t) := Xi{t) — x,(0) be the ith 
component of the displacement, we have 

(AxiAxj), ~ 2Dijt, (2.41) 

where the D,y are components of a symmetric diffusion tensor. 



(b) Central limit theorem: convergence to normal distribution Scale the displacement dis- 
tribution by so that the variance of the rescaled distribution is bounded. Then this distribution 
converges weakly, or in distribution, to a normally distributed random variable z IIGN901 IC YOOl : 

x(t)-x(0) V 

^ ^ z, as f ^ oo. (2.42) 
In the 1 -dimensional case, this means that 

^^LZ^<u]=—L= I t-^'l^^'ds, (2.43) 
y/t J aV27i J 

.V= — oo 

where P(-) denotes the probability of the event inside the parentheses with respect to the initial 
distribution of the random variable xq, and is the variance of the limiting normal distribution. 
In <i ^ 2 dimensions, this is replaced by similar statements about probabilities of J-dimensional 
sets. This is the central limit theorem for the random variable Ajc. From (a) we know that in ID, 
the variance of the limiting normal distribution is = 2D; md^2 dimensions, the covariance 
matrix of z is given by the matrix (2Z),y) HBSSIIIDCOOII . 
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(c) Functional central limit theorem: convergence of path distribution to Brownian motion 

Consider the following rescaling of the whole path of the process: 

„ , , x(st)-x(0) 

where s G [0, 1]. The scale y/t is the 'natural' scale coming from (b). This rescaling 'squashes' 
the entire path of the process onto the interval [0, 1], so that we can compare paths at different 
times. In fact, we compare the induced probability measures on the space of continuous functions 

[0, 1] ^ R'^. 

We say that the process satisfies afunctional central limit theorem if the probability distribution 
of the rescaled paths of the process converges weakly to Wiener measure (Brownian motion) B 
with covariance matrix as in (b), as f ^ oo; 

X; ^ B as f ^ oo, (2.45) 

This is known as afunctional central limit theorem, or weak invariance principle BCYOOII . 

This makes precise the sense in which X; looks like Brownian motion on long length and time 
scales. In the following sections we discuss in more detail the meaning of the above statement, by 
showing how diffusion can be regarded as a stochastic process and in what sense these rescaled 
processes converge to diffusion processes. 

2.4.1. Diffusion as a stochastic process 

As recalled in Sec. |1.3[ diffusion is described classically by the diffusion equation 

dp{t,r) 



dt 



:DV^p(?,r). (2.46) 



We would like a microscopic model which gives behaviour on a macroscopic level consistent with 
this equation. Following Einstein and Wiener (see e.g. ||Gar851 ). we look for a stochastic process 
Bt, determined by the probability density p{x,t) of a particle being at position x at time t given 
that it started at x = at time t = 0. Based on physical reasoning, we impose that the sample paths 
of the process should be continuous, and that displacements Bt+h — Bt should be independent of 
the history of the particle motion up to time t. 
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Under certain technical conditions, p{x,t) then satisfies the equation 

: 0, (2.47) 



dp d_ 

dt dxi 



2 ; dxj 

J •' 



which is known as Kolmogorov's forward equation or the Fokker— Planck equation UGarSSI . The 
drift vector A(x,f) and the diffusion tensor B{x,t) give the mean and variance, respectively, of 
infinitesimal displacements at position x and time t IIGar85ll . 

If the system is homogeneous, then A and B are independent of x and t. If the system is also 
sufficiently symmetric, then the drift is zero and the diffusion tensor is a multiple of the identity 
tensor. The stochastic process is then Brownian motion, and the Fokker-Planck equation (12.471 1 
reduces to the diffusion equation (12.461 ). (We remark that the factor 1 /2 in (12.471 ) occurs naturally 
in the probabilistic setting, and often leads to a discrepancy between the probability and physics 
literature.) A general diffusion process, however, can be inhomogeneous in both space and time. 

Brownian motion is defined as follows ||Dur96ll . A standard one-dimensional Brownian motion 
(also called the Wiener process) is a real-valued stochastic process Bf,t^O, such that: 

(i) For all n, if to < ti < . . . < t„, then the increments B{to),B{ti) — B(fo), • • • ,B{tn) — B(f„_i) 
are independent random variables. 

(ii) If t > s ^ 0, then B{t) — B{s) is a normal random variable with mean and variance t — s, 
so that 

/•I r r2 1 

dx. (2.48) 



V{B{t)-B{s) G A) = /■ ] exp 
J ^/2n t-s) 



A 



(iii) With probability 1, the function t ^ B{t) is continuous. 



l{t-s) 



Standard tZ-dimensional Brownian motion is then the vector process B(f) := (Si {t),.. . ,Bd{t)), 
where each B, is an independent standard ID Brownian motion. 

Another approach to study such diffusion processes is via stochastic differential equations 
(SDEs) IIGar85l . The above process corresponds to solutions of the SDE 

dx = A(x,f)df + a(x,OdB(0, (2.49) 
where B is an n-dimensional Brownian motion and <7(7^ = B. This form makes more explicit 
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the infinitesimal increments referred to above, and is now in a form which is more suitable for 
numerical simulation. 

2.4.2. Weak convergence to Brownian motion 

We must make precise what we mean by saying that 

Xt^B as f ^ oo. (2.50) 

The rigorous definition of this notion of weak convergence is given in App. see IIB 11681 . 
Necessary and sufficient conditions are [B 1168 II (1) that the finite-dimensional distributions of the 
process x, converge to those of Brownian motion; and (2) that the class of induced measures on 
path space is tight: see App. Ofor the definition. 

Property (1) means that for any n, any times < • • • < s,i, and any reasonable sets Di , . . . ,D„ 
in W', we have 

P (x,(5i) G Di, . . . M^n) G AO ^ IP (B(5i ) G Di , . . . ,B(5„) G D„) . (2.5 1) 

This is called the multi-dimensional central limit theorem in IIChe95ll . The right-hand side can be 
expressed as a multi-dimensional integral over Gaussians, as follows. We see thaj^ 

P (B(f2) = X2 I B(fi) = xi) = P (B(f2) - B(fi) = X2 - xi I B(fi) = xi) (2.52) 

= P(B(f2)-B(fi) = X2-xi) (2.53) 
= 7?(x2-xi,f2-fi), (2.54) 

since by the definition of Brownian motion B, the events B(f2) — B(fi) = X2 — xi and B(fi) = xi 
are independent, and B (?2 ) — B (f i ) is normal with mean and variance 2Dij {t2 — ti); here p (x, f ) is 
the Gaussian probability density function with that mean and variance. Integrating over all xi G D\ 
and X2 G D2, we have 

P(B(fi)GDi,B(f2)GD2) = j dxi j dx2/7(xi,fi);7(x2-xi,f2-fi), (2.55) 

Di D2 

^The probability that B(fi) = xi exactly is 0, but the argument can be made rigorous. 
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with a similar expression for the right hand side of (12.51b . 

Since x(s;) G D, if and only if x{tSi) — x(0) € \/tDi, the multidimensional central limit theorem 
becomes 

P(x(A2fi) G ADi,x(A^f2) G AD2,...,x(A\) G AD„) 

^ J dxi J dx2---J dx„p(xi,fi)p(x2-xi,f2-fi)---p(x„-x„_i,f„-f„_i), (2.56) 

Di D2 D„ 

setting X := \/t and then renaming the times si as This relation was used in MDCOOI as the 
definition of a diffusive process, although the functional central limit theorem is a stronger state- 
ment which also requires tightness. We remark that another name for the functional central limit 
theorem is the weak (or Donsker) invariance principle IICYOOI . 

In general, we may have convergence to a more general diffusion process than Brownian mo- 
tion. For example, in IIDGL81I it was shown that the motion of a large particle embedded in 
an infinite ideal gas converges to an Omstein-Uhlenbeck process. In the case of the periodic 
Lorentz gas, however, we have time-independent dynamics defined on a system which is space- 
homogeneous and symmetric on a large scale; the only possible limiting diffusion process which 
satisfies these conditions is Brownian motion. 

2.4.3. Further remarks 

An important part of the above limit theorems is proving that the limiting distribution is non- 
degenerate, i.e. that the diffusion tensor is positive definite IIBunOOII . In general, we can consider 
the above limit theorems for an arbitrary observable /. 

We will mostly consider 1 -dimensional projections (marginals) of the above distributions; we 
can equivalently regard this as studying the limit theorems for / being one component of Ax, 
which follow from the multi-dimensional results stated above. 

2.4.4. Discussion of definitions 

Property (c) makes precise in what sense a dynamical system looks like Brownian motion when 
correctly rescaled. This is the strongest, and so in some sense the best, property with which we 
could define deterministic diffusion (i.e. a dynamical system is diffusive when it satisfies property 
(c)). However, there are very few physically relevant systems which have been proved to satisfy 
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(c). Interest in the periodic Lorentz gas comes in large part from the fact that it is one of the 
only such systems; another is the triple linkage IIHM03II . As mentioned above, in general we may 
instead have convergence to some other diffusion process. 

The multi-dimensional central limit theorem part of (c) was studied in MDCOOi . where both 
Lorentz gases and wind-tree models were found to obey it, tested for certain sets D, and certain 
values of n. However, as stated in BDCOOI . (c) is difficult to investigate numerically, and the results 
in that paper seem to be the best that we can expect. 

Property (b), the central limit theorem, has been shown for large classes of observables / in 
many dynamical systems (see e.g. IIChe95ll and references therein), but again they are not often 
physically relevant systems. Property (b) was used in IIGN90II as the definition of a diffusive 
system, but does not seem to have been applied in the physics literature; this is the approach taken 
in Chap, m 

Most emphasis in the physical Uterature is placed on property (a): many authors define a system 
to be diffusive if only property (a) is verified (numerically), e.g. IIKDOOi lA"RdV02[ iDCOl I . Many 
types of system are diffusive in this sense, including ID maps IIKD99II . Lorentz gases (periodic 
and random) MBSSlilKDOOllDCOll and Ehrenfest wind-tree models, both periodic IIARdV02ll and 
random liDCOll . 

The implication (c) =^ (b) always holds, and (a) and (b) usually go together. The reverse im- 
plications (a) =^ (b) and (b) =^ (c) hold only under certain additional conditions: see e.g. IIChe95l . 
It is thus of interest to ask if dynamical systems showing property (a) also show properties (b) and 
(c). For example, in IIKC89i a disordered lattice-gas wind-tree model was reported as having an 
asymptotically linear mean squared displacement, but a non-Gaussian distribution function, i.e. (a) 
but not (b). However, disorder can lead to trapping effects which cannot occur in periodic systems 
IIARdV02ll ; we are not aware of a periodic (and hence ordered) billiard model with a unit-speed 
velocity distribution which shows (a) but not (b), although in Sec. 14.71 we show that this can occur 
with a Maxwellian velocity distribution. 

2.4.5. Approach via cumulants 

Several papers have approached property (b) by studying higher order cumulants of the distribu- 
tion, e.g. IIARdV021IDC01ll . A normal distribution has all cumulants higher than the second equal 
to zero. It thus seems that we need all cumulants of the rescaled displacement distribution to tend 
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to zero as f ^ oo for the distribution to approach a normal distribution; this was stated in HDCOU 
p. 790], for example. 

In fact, the weak type of convergence occurring in the central limit theorem (convergence in 
distribution) does not in general require higher moments and cumulants even to exist: a faster rate 
of decay of correlations is necessary to ensure the existence of higher cumulants than is required 
for the central limit theorem to hold: see llvB82[|CD00l . 

Higher cumulants are instead related to large deviations, describing behaviour in the tails of 
the distribution IIGas98[ Sec. 7.3.6]. We show in Chap.[5]that a class of polygonal billiard channels 
which was found in IIARdV02ll to have rapidly-growing cumulants nonetheless appears to satisfy 
the central limit theorem. 

2.5. Rigorous results on billiards 

In general, a billiard table is a Riemannian manifold with a piecewise smooth boundary. The 
dynamics is given by geodesic flow away from the boundary and specular reflection (elastic colli- 
sion) at the boundary. Here we restrict attention to two-dimensional periodic Lorentz gases, where 
non-interacting point particles in free motion on a torus-shaped billiard table undergo elastic col- 
lisions with a set of fixed scatterers. We give only enough detail for our purposes, and refer to 
IITab95L liCFS821 Chap. 6] and liCMOli Chap. 4] for further informatioif . 

2.5.1. Continuous-time dynamics 

We denote by Q the billiard domain, i.e. the available region where particles can move, given by 



where the scatterers Vi are non-intersecting discs in the case of the periodic Lorentz gas. We 
visualise the torus as a square [0, 1)^ with periodic boundary conditions, as in Fig. 12.21 

Each particle moves in a straight line within Q at constant velocity v until it hits the boundary 
dQata point q on the boundary dVi of scatterer /. It then undergoes an elastic collision, giving a 
new velocity 



(2.57) 



v' = V- 2(v-n)n 



(2.58) 



^[An excellent recent reference on chaotic billiards is ICM06I .1 
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Figure 2.2: Example of a billiard domain Q with scatterers removed. 



after collision, where n := n(q) is the outward normal vector at q to the scatterer P,. The particle 
continues its motion by alternation of straight line motion and elastic collisions with the boundary. 

The above prescription defines a continuous-time dynamical system given by a flow 

^'■.M^M; (2.59) 
a>^(qo,Vo)^(q„v,), (2.60) 

where A4 := TQ is the tangent bundle to Q, i.e. the set 

M:={{q,y):q£Q,y£%,Q}, (2.61) 

where T^Q is the tangent space to Q at the point q. 

Since the speed ||v|| is preserved, it is customary to take all particles with unit speecj^ We can 
thus restrict attention to the new phase space M given by the unit tangent bundle ^ 

M:=T^Q = QxS\ (2.62) 



where 

5^ := {ve : ||v|| = l}, (2.63) 

and the flow ^ : M Ai, which we call the billiard flow. We can parametrise the space Ai with 
coordinates {x,y, 6), where {x,y) G [0, 1)^ are position coordinates and 6 G [0,27r) is the angle of 

^We consider non-constant speed distributions in Sec. 14. 71 

^The identification of the unit tangent bundle Ti Q with the direct product Q x 5' is not in general valid, but it does 
hold for the simple situation we consider. 
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the velocity vector. 

We remark that for billiard tables with comers, where two smooth parts of the boundary touch 
at a non-zero angle, the trajectory cannot be continued if it exactly hits a comer. However, the set 
of all such trajectories has measure 0. 

2.5.2. Billiard map: Poincare section 

As is often the case in dynamical systems, it turns out to be technically easier to study a Poincare 
section of the flow. A natural cross-section is given by 



i.e. the set of vectors whose base points lie on a scatterer and which point into the interior of the 
billiard domain. This cross-section M forms a two-dimensional surface in the three-dimensional 
phase space Ai. One way of parametrising M is to take as coordinates the arc length r around 
the scatterer of the point of collision (from some reference point), and the angle (p between the 
velocity vector and the normal vector IICMOlll . Altematively we can take the arc length around the 
scatterer and the sine of cp, giving canonically conjugate Birkhoff coordinates liGas98[ Chap. 5]. 
In general we must also keep track of a label / denoting which scatterer was hit. The flow then 
induces an invertible map T : M ^ M which maps one collision into the next, called the billiard 
map. 

The billiard flow can now be regarded as a suspension flow over T, under the free path 
length function x :M ^ M+, where T(q,v) is the length of the collision-free part of the trajectory 
emanating from the initial condition (q, v) G M. Suspension flows are treated in detail in App.lBl 

2.5.3. Measures 

Since the biUiard flow is Hamiltonian, it preserves a natural invariant measure }X, called Liou- 
ville measure, given 



M ■={{q,\) £ M: q e ^2, v • n(q) > 0} 



(2.64) 



d/i := c^dxdydd, 



(2.66) 



This is a useful shorthand for the statement 




for all AeB. 



(2.65) 



A 



(x,y,e)eA 
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where dx etc. denote Lebesgue measure in the respective coordinate, and is a normalising 
constant chosen so that IJ.{M.) = 1. This measure induces a measure v on M which is invariant 
under T, given by MCMOll 

dv := CvCos(pdrd<p. (2.67) 



For a 2D billiard (i.e. one for which the dynamics is restricted to a plane), we have 

1 



^ 2n\Q\ 

and 



(2.68) 



2.5.4. Ergodic and statistical properties 

The ergodic and statistical properties of billiards have been the object of intense scrutiny since the 
seminal work of Sinai liSinTOl . where the K-property of the periodic 2-disc fluid was proved. All 
of the proofs depend on the hyperbolicity of scattering billiards, i.e. the existence of stable and 
unstable directions; for details see [ CM06[ ICMO 1 1 lGas98[ ITab95i and the references cited below. 
A few of the key rigorous results for dispersing and semi-dispersing billiards are: 



• 2D periodic Lorentz gas models for which a geometrical finite horizon condition holds 
(Sec. 13.1.11 ) satisfy the central limit theorem and functional central limit theorem if the 
scatterers are disjoint and piecewise smooth IIBS81[|BSC91I : 

• higher-dimensional Lorentz gases with the same conditions also satisfy the CLT and FCLT 
IIChe941 . but there are issues with the proofs in higher dimensions that need addressing^ 
after the work of IIBCST03[|BCST02II : 

• velocity autocorrelation functions of the billiard map decay exponentially IIYou981IChe99ll : 
and 

• hard ball fluids are ergodic, mixing and K-sy stems: see llSimI and references therein. 

There is an excellent collection of reviews in llSzaOOII . 
"'N. I. Chernov, private communication. 
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2.6. Numerical evaluation of statistical quantities 

To evaluate numerically statistical quantities such as the mean squared displacement, we use a 
simple Monte Carlo method. We take a large sample (xq^Vq^)^! of size of initial conditions 
chosen uniformly with respect to Liouville measure in one unit cell using a random number gen- 
erator: the positions xq are uniform with respect to Lebesgue measure in the billiard domain Q, 
and the velocities vq are uniform in the unit circle S\ i.e. with angles between and 2n, and unit 
speeds. These evolve after time t to (x^'^ (t), y^'^ (0 > the distribution of this ensemble then gives 
an approximation to that of {x{t),\{t)). 

We use the Mersenne Twister random number generator IIMN98I . which seems to be consid- 
ered to be one of the best-performing generators; we have also done some comparisons with a 
standard rand implementation. Since our random numbers determine only the initial conditions, 
and since the system is strongly mixing, we do not expect the choice of random number generator 
to be important, provided the initial conditions are sufficiently well distributed through the phase 
space, so that no small regions with important dynamical effects are omitted. 

We denote averages over the initial conditions, or equivalently expectations with respect to 
the distribution of (xo,Vo), by (•). Approximations of such averages can now be evaluated via 
IIPTVF92I 




(2.70) 



The infinite sample size Mmit, although unobtainable in practice, reflects the expectation that larger 
N will give a better approximation. 

Lyapunov instability The Lyapunov instability in dispersing billiards implies that numerical 
simulations of trajectories will not be accurate beyond a time where the uncertainty in the initial 
condition has grown exponentially to 0(1), so that some authors do not allow any data obtained 
beyond that time IIGG94II . However, most physics papers implicitly reject this argument by pre- 
senting long-time results from simulations, and we shall also do this, whilst being careful about 
other sources of error in statistical calculations. One partial justification for this comes from shad- 
owing IICM06I : in certain hyperbolic systems, it can be proved that there is a true trajectory close 
to a simulated trajectory. The difficulty in billiards comes from the discontinuous nature of the 
dynamics; nonetheless, some results of shadowing type were proved in IIKT92II . We are not aware 
of any papers which investigate shadowing in particular billiard models, but early numerical in- 
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2. Statistical properties of dynamical systems 



vestigations of the validity of numerical estimates of ergodic propeities for uniformly hyperbolic 
systems can be found in IIBCG+TSIIbCG+TQII . 



Chapter 3 



Geometry-dependence of the diffusion coefficient in a 2D 
periodic Lorentz gas model 



3.1. Two-dimensional periodic Lorentz gas model 

A periodic Lorentz gas is a periodic billiard with disjoint circular scatterers (discs). The term is 
also applied to other periodic billiards of dispersing type with scatterers of more general shape: 
see e.g. IIHKG02I . 

In this chapter and the following one we study aspects of deterministic diffusion in a partic- 
ular 2D periodic Lorentz gas model introduced below: here we consider only the mean squared 
displacement and the diffusion coefficient, while in Chap. |4] we examine the shape of distributions. 

3.1.1. Finite horizon condition 

Periodic Lorentz gases were shown in IIBS811|BSC91|| to be diffusive (Sec. l2.4l ). provided they sat- 
isfy the finite horizon condition: there is an upper bound on the free path length between collisions. 
If this is not the case, so that a particle can travel infinitely far without colliding with any scatterers 
(the billiard has an infinite horizon), then corridors exist, which allow for fast-propagating orbits. 
This leads to super-diffusive behaviour, in the sense that the mean squared displacement grows like 
tlogt (which is faster than t), as has recently been proved in IISV07I . after heuristic and numerical 
arguments IIFM841 IZGNR861 and analytical evidence IIBle92n were previously given. There is a 
more detailed discussion of statistical behaviour in the infinite-horizon regime in Sec. 16.41 

We aie mainly interested in normal diffusion, so that we wish to restrict attention to values 
of the geometricalnarameters for which there is a finite horizon. The simplest periodic two- 
dimensional latticqj is the square lattice. To allow the possibility of diffusion in a square lattice, 



'We consistently use the term 'lattice' in the sense of physics to mean any periodic structure. Mathematicians 
instead use the term in the narrower sense of the set of points in of the form '^i^i^ where the e,- are basis vectors 
of and a; £ Z; these are called Bravais lattices in physics IAM76I . 
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3. Geometry-dependence of diffusion coefficint 




Figure 3.1: (a) Part of the infinite system, constructed from two square lattices of discs, shown in different 
shades of grey; the dashed lines indicate several unit cells, (b) A single unit cell, defining the geometrical 
parameters. The billiard domain is the area Q exterior to the discs. 

the scatterers must not touch or overlap: if they do, we obtain a localised regime where any particle 
is trapped in a bounded region of space, so that the diffusion coefficient is necessarily 0. However, 
non-touching scatterers in a square lattice always imply an infinite horizon: trajectories parallel to 
the lattice axes and passing sufficiently close to the centre of a unit cell never collide. 

For this reason a periodic Lorentz gas on a triangular lattice is often used, where there is a 
finite horizon when the non-touching discs are sufficiently close together: see e.g. IIMZ83[|KD00I . 
We instead elect to retain a square lattice model, blocking all possible corridors by adding an 
additional disc at the centre of each cell. 



3.1.2. Model studied 

The model which we focus on, previously studied in ||GG941lGar97l . consists of two inter-penetrating 
square lattices of discs; they have the same lattice spacing r, and radii a and b, respectively, and 
are arranged such that each unit cell of the a-lattice contains a Z?-disc at its center. Part of the 



infinite system is shown in Fig. 3.1(a) and the geometrical parameters are defined in Fig. 3.1(b) 



We were led to this model as a cross-section of the 3D model presented in Chap. [6l indepen- 
dently of ||GG941lGar97ll . although as described above, it is a natural model to consider. The same 
model has also been studied in IIBDLR021 IBDLOOII in a different context. The two-dimensional 
parameter space of this model allows us to vary independently, and hence study the effect of, two 
physical quantities: the size of exits of each unit cell, and the accessible area of a unit cell. This is 
not possible in the standard triangular Lorentz gas. 



3.1. Two-dimensional periodic Lorentz gas model 
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3.1.3. Length scale 

Since the speed is conserved, scaling the geometry by a constant factor does not alter the essentials 
of the dynamics, just the time-scale on which it occurs. We may hence use dimensionless units 
such that the speed is equal to 1 and we are free to choose a length scale. 

In this thesis we use one of two choices of length scale: in analytical calculations it is often 
convenient to fix a = 1 and vary r and b; while in numerical calculations we often fix r = 1 and 
vary a and b. Where necessary we distinguish the latter case using tildes, writing 

d:=a/r, b:=b/r. (3.1) 



3.1.4. Parameter space 

Garrido IIGar97l derived the 'phase diagram' of the model, showing the regions in parameter space 
corresponding to the various localised, finite horizon and infinite-horizon regimes. Here we correct 
this diagram by including the following features that were incorrectly treated in ||Gar97ll : (i) the 
system is symmetrical under interchange of d and b; and (ii) the finite-horizon regime is smaller 
than was found in ||Gar971 . We fix r = 1 but omit the tildes on d and b. The derivation, which 
reduces to geometrical considerations, is useful since it provides pointers for the geometrically 
more involved calculation for the 3D model in Chap. |6l 

If we interchange the radii a and b, then the resulting structure is the same as the initial one, 
except for a translation by the vector {\,\)- Hence the statistical properties are the same, so 
that the whole phase diagram is reflection symmetric in the line a = b.\n the following we hence 



restrict attention to the triangular region a'^b shown in Fig. |3.2(a)| and then complete the diagram 
by reflection. 

The main features of the geometry are (i) whether the a discs touch or overlap with each other; 
and (ii) whether the b discs overlap with the a discs. The respective boundaries in parameter space 
are given by the straight lines with equations 

£2 = 1/2 and b = -^ — a, (3.2) 

' \/2 

respectively. These divide the region a> b into 4 sub-regions A, B, D {diamond) and S {star), 
shown in Fig. |3.2(a)| 
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In the complete phase diagram, shown in Fig. |3.2(b)| we have the following subdivisions: 

• region A into FH (finite horizon) and IH (infinite horizon); 

• region B into T (triangle) and O (overlapping: whole plane covered). 

The IH regime can be further subdivided into regimes with different classes of possible infinite 
trajectories: 

• IHl: both vertical/horizontal and diagonal infinite trajectories exist; 

• IH2: only diagonal ones exist; and 

• IH3: only vertical/horizontal ones. 

The region IHl can itself be subdivided into infinitely many regions where increasingly many 
corridors (at diagonal angles other than 45°) become available. Bleher ||Ble921 Sec. 8] gives a 
concise method to find the boundaries of these different IH regimes in the case b = 0. 

To have a finite horizon we must block horizontal, vertical and diagonal corridors. (The neces- 
sity to block diagonal corridors was overlooked in ||Gar97l| .) Horizontal and vertical corridors are 
blocked by requiring the central disc to be sufficiently large, namely b ^ ^ — a, and by symmetry 
it suffices to block diagonal trajectories at 45° between mid-points of adjacent sides; for this the 
a discs must be large enough. (If we instead use b discs to block these diagonal trajectories, we 
necessarily have b>a;by our convention, we thus reduce to the case considered.) Hence for a 
finite horizon we supplement Garrido's condition 

^-a ^ b ^-a (3.3) 

witlfl 

a>-^. (3.4) 

The boundary between regions T and O occurs when the b disc covers the whole region left 
empty between the overlapping a discs. This occurs when the boundary of a ft disc reaches the 
point of intersection of the boundaries of two a discs, when ft = ^ — w — l. 



^Note that the critical case of a 'grazing' collision (when a trajectory is tangent to a scatterer) is sufficient to have 
normal diffusion. 



3.1. Two-dimensional periodic Lorentz gas model 
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Figure 3.5: One-dimensional channel created by unfolding torus only in jc-direction. 



The complete phase diagram is shown in Fig. 3.2(b) The dots give the parameter values used 



to depict the various regimes in Fig. 13.31 (non-localised configurations) and Fig. 13.41 (locahsed 
configurations); white indicates the allowed region of motion of the particles. 

3.1.5. Channel geometry 

Diffusion in a 2D infinite lattice is described by a second order diffusion tensor having 4 compo- 
nents Dij with respect to a given orthonormal basis, which we refer to as diffusion coefficients. If 
the system is sufficiently symmetric then the diffusion tensor reduces to a scalar multiple of the 
identity tensor. Square symmetry, as in the model presented above, and hexagonal symmetry, as in 
the standard triangular periodic Lorentz gas IKDOOl are both sufficient for this reduction to occur; 
see also Sec.f 



We study individually the components of the diffusion tensor, defined by 

Dij = lim ^ {AxiAxj)r . (3.5) 

This enables us to check that D^^ = Dyy =: D and D^y = Dy^ = in fully symmetric systems. For 
those systems the diffusion coefficient can then be evaluated as an average over the multidimen- 
sional distribution via 

(Ax^), = (Ax^), + (A/), ~ ADt, (3.6) 

as used e.g. in MKDOOI . but this cannot be applied to systems where the diffusion tensor is not a 
multiple of the identity tensor. 

We can evaluate the diffusion coefficient D^x by looking at the dynamics only in the jc-direction, 
which corresponds to studying the billiard dynamics in a 1 -dimensional channel extended in the x- 
direction; see Fig. 13.51 Correspondingly, we restrict attention to ID marginal distributions, which 
are easier to analyse. 

A channel geometry, with hard horizontal boundaries, corresponding to the triangular Lorentz 
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(b) 



Figure 3.6: (a) Lorentz channel studied in MGas93l |AACG99]| with hard upper and lower boundaries; 
dotted lines indicate unit cells, (b) Fully unfolded triangular Lorentz gas. Dotted lines indicate unit cells 
forming a channel with periodic upper and lower boundaries, created by reflecting the channel in (a) in the 
upper boundary. 



gas was studied in IIGas931 lAACG99i (Fig. |3.6(a)[ ). This is equivalent to a channel with twice the 
original height and periodic boundaries, by reflecting once in the hard boundary. This new channel 
is shown in Fig. |3. 6(b)] as part of the whole triangular lattice obtained by unfolding completely in 
the vertical direction. We can then view the triangular lattice as consisting of rectangular unit cells 



(Fig. |3.6(b)| ) which are stretched versions of the square unit cell considered above, with the extra 
condition a = b. 



3.1.6. Effect of length scale on diffusion coefficient 

Since we have two choices of length scale (Sec. |3.1.3l l. we need to know how to convert calculated 
diffusion coefficients between the two length scales. Recall the convention that d = a/r, b = b/r; 
we also introduce for clarity r := r/r = 1. 

Let x{t) be the position at time t in the system with a = I fixed and speed 1, with radii r and b; 
and let x{t) be the position at time t in the system with f = 1 fixed and speed 1, with radii d and b. 

The motion in the tilde system is then equivalent to that in the a = I system, but with speed r; 
however, the size of the system is also shrunk by a factor r. Thus we have 




(3.7) 



and so 



<*« = (C-y^)'')~>'^ 



(3.8) 
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Hence 




(3.9) 



where D is the diffusion coefficient in the tilde system. (This result can also be obtained by 
dimensional analysis.) 

3.2. Estimation of diffusion coefficient 

Having set up the model, we now consider the statistical question of how to estimate the diffusion 
coefficient from numerical data. 

3.2.1. Estimation of moments 

A first approach to characterise the distribution of a random variable is to look at its moments. 
The rth raw population moment of a random variable X is given bjcl 



where Fx is the distribution function of the random variable X ; the integral with respect to Fx is a 
Lebesgue-Stieltjes integral, and the last equality holds if the random variable has a density fx- 

If we take a sample (X,),=i a? of independent and identically distributed random variables 
from the distribution of X, then the mth sample raw moment is 



This is an unbiased estimator of the rth population raw moment jx'^, i.e. we have E [m'^ = jx'j.. 
Halmos IIHal46ll showed that in fact it is the unique unbiased estimator of jU^(X) which is a sym- 
metric function of the X,, and further that this estimator has the smallest variance of all unbiased 
estimators; in this sense it is the best estimator of IJ.^{X) given the sample. 

We are mainly interested in the distribution of position and displacement in billiard models 
starting from an initial distribution which is uniform in a unit cell. By symmetry, the mean dis- 
placement {Ax)t then always vanishes, so that the simplest non-trivial characteristic of the dis- 

%e follow the notation of IRS02I . 




(3.10) 



/ 1 

m^ := — 

' N 



N 



(3.11) 
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tribution is the mean squared displacemen\\ (Ax^),. We can view (Ax^), either as the 2nd raw 
moment of the random variable Ajcj, or as the 1st moment of the distribution of the new random 
variable Yt := Ax^^, as follows: 

{M%=^'^{^Xt)=^[{Yt), (3.12) 
The best estimator of the mean squared displacement is thus 

m\ := m\{Y,) = m^(Ax,) = i (3-13) 

this can be regarded as a simple Monte Carlo estimator (see Sec. 12.61) . In our numerical exper- 
iments we restrict ourselves to a fixed number of initial conditions N, although in principle we 
could add more data until some pre-defined error tolerance was reached. 



3.2.2. Distribution of mean squared displacement at time t 

To estabhsh how good an estimator (13.131 ) is, we need to determine the width of the distribution 
of m\. The ^sxf' are independent and identically distributed (i.i.d.) random variables with mean 
0; since particles with speed v satisfy \^t\ ^ vt, they also have finite variance. It follows that the 
7(0 := (Ax|'^)^ are also i.i.d., with positive mean and finite variance. Hence the standard central 
limit theorem (Sec. 12.3.21 ) applies, so that for large N the distribution of m\ is very close to normal, 
with unknown mean and variance 

Var[m;] =i^IVar[Zk?] =ivar[Zk2]. (3.14) 

For fixed t, we thus arrive at the standard result that the width of the distribution of m'j , as measured 
by the standard deviation (the square root of the variance) of m\, is 0{N^^^^). 

The best estimator for this variance can be found using the known result for the variance of a 
sample mean, giving 

S':=j^^t(Y}''>-m[iY,))\ (3.15) 



^The notation (Ax^), emphasises that we are averaging over the distribution at time t, but we have {Ax^)i = 
{[Ax{t)]'^), where now we are thinking of averaging over the distribution of initial conditions. This distinction is the 
same as that between the Schrodinger and Heisenberg pictures in quantum mechanics. 
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3.2.3. Calculation of error bars: confidence interval for {^^^)t 

We now establish a confidence interval for (Ax^);. For a normal random variable Z with mean pi 
and variance a we have 

P e [2.576,2.576]^ = 0.99, (3.16) 

i.e. the probability that a normal random variable lies within 2.576 standard deviations of the mean 
is 99%. Since neither }JL nor a are known for the distribution of ^[iYt), we estimate them via m\ 
and s, respectively, where s is the value of for the data. 

We will then estimate a 99% confidence interval for fi'^Yt) by the interval estimator 

[m'l - 2.5765, m'l + 2.576^] . (3. 17) 

Figure [377] shows the estimate m\{Yf) as a function of t for several samples. Error bars for one 
sample are also shown at m\ it 2.576^; the error bars contain most of the data for each sample. 

3.2.4. Estimation of diffusion coefficient 

The most efficient method to estimate the diffusion coefficient is to use the fact that the asymp- 
totic growth rate of the mean squared displacement in the linear regime is 2D and to calculate 
this by fitting a straight line to the data in the presumed linear regime IIKDOOI . We must thus 
consider the evolution of (Ax^); over time. In practice, the linear regime is attained very rapidly: 
see Fig. I3.7f a). To find the slope in the asymptotic regime we apply linear regression (see e.g. 
IIPTVF92I ). performing a least squares fit of a straight line to the data in the region of large t, 
making sure that the correlation coefficient of the fit is very close to 1 . 

The Unear regression procedure gives an estimate of the diffusion coefficient as the slope of 
the fitted line. We also need to estimate the possible error that this estimate has made compared 
to the true value. This question is treated in statistics textbooks mostly for data whose errors are 
independent. This is however not the case here: if at a given time t the estimated mean squared 
displacement has deviated from the 'correct' value, then at a small time later, there is a large 
probability that the mean squared displacement deviates away from the 'correct' value in the same 
direction as in the previous step. 
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Figure 3.7: (a) Mean squared displacement (Ax^), as a function of t for 20 samples, each averaged over 
A' = 10000 data points. Also shown are error bars showing 99% confidence intervals for one of the samples. 
The other data mostly lie inside the error bars, (b) Deviations from a linear fit to the data set whose error 
bars are shown. 
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3.2.5. Growth of width of distribution of hx} with t 

To find the growth rate of the width of the Ax^ distribution, we argue as follows. Since the distri- 
bution of docf at a fixed time t is (approximately) normal, we only need to consider the standard 
deviation. We have 

Var [Axl] = E [Axf] - (E [/^] f = n'^{t) - n^{tf, (3.18) 

setting At,;(0 :=At,UA^f)- 

Chernov & Dettmann MCDOOH showed that for a periodic Lorentz gas with finite horizon, the 
4th-order Burnett coefficient B exists, where 

B:=\xm^K^{t), (3.19) 

and 

K4(0:=At4(0-3At2(0' (3-20) 



is the 4th-order cumulan^ This cumulant is related to the non-dimensional kurtosis excess 72(f) 
via 

the kurtosis excess 72(f) is a common measure of how far a distribution is from a Gaussian, the 
difference being exactly for a Gaussian liWeiH . 

We thus have 

yai[Axj] = K4{t) + 2n2{tf ^t + t^ r^t^ asf^oo, (3.22) 

since jJ.jit)'^ grows like whereas K4(f) grows like t. Hence the width of the distribution of Axf, 
as measured by the standard deviation, grows like t. This does not seem to have been previously 
remarked; it can be seen in Fig. [331 where the half- width 2.576 ^(Ajc^) of a 99% confidence 
interval for {Ax'^)t is plotted as a function of t. 

We note that Garrido & Gallavotti IIGG94II considered the question of estimating the error 
in the diffusion coefficient by an expansion assuming that the errors are small. However, if we 



^The set of cumulants of a distribution is an alternative set of moment-type descriptors of a distribution which have 
certain advantages over the raw moments: see e.g. IRS 021 . For example, K„{X + Y) = K„{X) + K„{Y) for independent 
random variables X and Y; this does not hold in general for mj,. 
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Figure 3.8: Growth of half-width of 99% confidence interval for (Ajc^), for two different sample sizes A^. 

could continue the simulation indefinitely, the above calculation shows that the errors would grow 
arbitrarily large. Nonetheless we can still estimate an error in the diffusion coefficient, as follows. 

3.2.6. Sampling distribution of diffusion coefficient 

By repeating the above estimation procedure for M different sets of random initial conditions, we 
obtain a sample set {D^'^)^i of estimated values of the diffusion coefficient. From these we can 
estimate the sampling distribution of D; this will indicate how close to the true underlying value 
of D the estimate D^'^ is expected to be. 

Suppose that the confidence intervals for the mean squared displacement enclose all of the 
data. Then the growth rate of the data must lie between the rates of growth of the upper and 
lower limits of the error bars. In fact we can never guarantee that all the data is enclosed, but 
nonetheless we estimate the width of the D distribution by fitting straight lines to the error bars. 
For 99% confidence intervals on the mean squared displacement error bars, we would naively hope 
to obtain a 99% confidence interval on D. However, a confidence interval for each point does not 
necessarily correspond to a confidence band for the whole curve ||Ric94|| . 

We investigate the sampling distribution of D by constructing a kernel density estimate IISil86l 
of the sampling density. Fig. 13.91 shows the samphng density obtained for particular geometrical 
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Figure 3.9: Sampling density of D for M = 1000 samples of size = 10000, r = 2.5,b = 0.5, obtained by 
a kernel density estimate. A Gaussian with the same mean and variance is shown for comparison. 

parameters for M = 1000 samples of size N = 10000. For this case, the standard deviation of the 
sampling distribution is a = 0.044. The supposed 99% confidence interval found by the above 
method (fitting upper and lower error bar limits) is 0.08 ±0.012. (This value varies depending 
on the sample.) The estimated error bars from D hence correspond to a width of 1.8a, which is 
approximately a 90% confidence interval. Our method thus slightly underestimates the width of 
the sampling distribution. We have repeated the calculation for larger values of N, and have found 
empirically that our method of estimating the width of the distribution improves as N increases. 

3.3. Geometry dependence of diffusion coefficients 

We now consider the geometry-dependence of the diffusion coefficient in our model as a function 
of the geometrical parameters. This question has been studied in BKDOOI for the triangular Lorentz 
gas, and in IIHKG02II in a billiard with 'flower-shaped' scatterers. As emphasised above, our 
model has the advantage that we can vary independently two physical factors which we expect 
to influence the diffusion coefficient: (i) the size of trap exits; and (ii) the available area in each 
unit cell. We remark that striking results have also been obtained for lifted circle maps on the real 
line IIKD951IGK02II . where a fractal parameter dependence of the diffusion coefficient was found. 
This motivated a conjecture of Klages MKDOOII that the diffusion coefficient in low-dimensional 



3.3. Geometry dependence of diffusion coefficients 57 

dynamical systems could in general be a highly non-trivial function of parameteo We have not 
attempted to address this question: although we find that the diffusion coefficient in our model is a 
reasonably smooth function of parameter, we have not studied the fine-scale dependence in detail. 



3.3.1. Constant trap exit size 

In this section we work in the tilde scaling, i.e. fixing the side of the unit cell to be of length 1. 
The parameters varied are then a and b. 

We fix the radius a of the non-central discs and vary the central disc radius b over the allowed 
portion of the finite horizon regime for that value of a. We repeat this for different values of d 
covering the finite-horizon regime. Each curve thus corresponds to fixing the size of the exits of 
the unit cell, while the trap area varies. Figure ITTOl' a) shows the diffusion coefficient a) over 
the whole finite-horizon regime. The curves are plotted over the allowed range 

1 ~ 1 

--a^b<-j=-a=:bmax{a); (3.23) 

D vanishes for S ^ 1/ \/2 — d, since the particle is then completely trapped in a bounded region, 
so the rightmost symbol on each curve is plotted at (^max(^)>0). Error bars, estimated as in the 
previous section, are plotted, and are indistinguishable from the data symbols (except possibly for 
the largest values of d). 

We note the following features of these graphs. 

(i) The curves seem to be continuous as a function of b, with some degree of regularity. 

(ii) They possess features such as local minima and maxima which are reproduced in nearby 
curves, indicating that there is also continuity of D as a function of d with b fixed (see also 
below). 

(iii) Decreasing d with b fixed results in faster diffusion. 

(iv) There is a critical value a = — 0.48, below which the curves acquire a non-trivial maxi- 
mum away from the left end of the curve. Figure I3.10f b) shows in detail the region where 
this transition occurs, i.e. the lower left portion of Fig. [3?T0l 

^[ Rigorous results on this question have recently been obtained for simple maps IKHK08I and Lorentz gases ICDI 
chapter 5].] 
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Figure 3.10: (a) Geometry dependence of diffusion coefficients over whole finite-horizon regime. Each 
curve is for a fixed value of 5 = l/r, where 2.05 ^ r ^ 2.8 and r changes in increments of 0.05; r = 2.01 is 
also shown, (b) Detailed view of the region 2.01 ^ r ^ 2.1 in steps of 0.01. In both (a) and (b) r increases 
from bottom to top. 
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(v) For a close to its upper bound 1 /2, i.e. when the a-discs are very close together, the maxi- 
mum diffusion rate occurs for a value of b close to Z^max- 

Point (lulb corresponds to the findings with the triangular periodic Lorentz gas ||KK02i . since 
there both trap area and trap exit size increase when the parameter is increased. Intuitively, it is 
'easier' for the particle to move through the lattice. 



Point ^3 is reminiscent of the behaviour near to a critical point in the context of phase transi- 
tions. We have not found a physical explanation of this, although Sec. 13.41 contains some related 
comments. 



Diffusion coefficient as function of gap size Since the finite-horizon regime is a parallelogram. 



it is interesting to consider the diffusion coefficient in terms of a new variable c := S — (i — a) 



giving the distance above the least allowed value for b. In fact, we also have 

c = {^-\)-W2. (3.24) 

where W2 is the minimum distance from the boundary of the b-Ahc to the boundary of an a-disc 
(see also the next section). This is shown in 13.111 We can thus view the graph as plotting the 
diffusion coefficient as a function of this distance. 

We see that the graphs approximately collapse for values of W2 close to its maximum, i.e. when 
the a- and Z?-discs are close to touching. In particular, there seems to be a universal approach of 
the diffusion coefficient to when W2 0. We remark that it was shown in IIBun85ll that viewing 
D as a function of W2, we have 

^ const, (3.25) 

W2 

SO that this approach occurs in a linear way as D ~ W2; this should be contrasted, for example, 
with the square root- type behaviour D ~ ^/w2 often found in phase transitions. 

3.3.2. Diffusion coefficient variation with constant trap area 

We now fix the available area per unit cell and vary the trap exit size; again we work in the tilde 
system. Figure 13.121 shows the variation of the diffusion coefficient along circular curves in the 



parameter space of Fig. |3.2(b)| with constant radius a = \l cP- + b'^, and hence constant billiard 
domain area |2| = 1 — na^. The diffusion coefficient is plotted as a function of d := tan^'(^/a), 
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0.14 




Figure 3.11: Geometry dependence of diffusion coefficients over tlie whole finite-liorizon regime, as a 
function of c :— b ~ {j ~ a). Tfie curves are as in Fig. 13.101 

the angle anti-clockwise from the ^-axis along the circle of radius a. When we move along these 
circles the radius b of the central disc increases, while the radius d of the other discs decreases, 
and hence so does the trap exit size. We find the following: 

(i) As a decreases, and so the available area increases, the rate of diffusion increases. 

(ii) For larger values of a, the curves lie mainly within two regions where a random walk model 
gives a good approximation to the diffusion coefficient, as we show in Sec. 13.41 

(iii) For a « 0.45, the main effect seems to be due to variation of the trap area. However, there is 
a still some variation along the curves. This provides some evidence against the application 
of an approximation originating in the Boltzmann equation proposed in HKDOOL since that 
approximation depends on density alone. In that paper the application was to the triangular 
periodic Lorentz gas, but the lack of variable parameters in that model did not allow the 
investigation of this question. 

3.4. Machta-Zwanzig random walk model 

We wish to understand the dependence of the diffusion coefficient D on the geometry of the pe- 
riodic Lorentz gas: can we predict the gross structure of curves such as those in Sec. 13.31 via 
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Figure 3.12: Diffusion coefficient as a function of angle 9 around circles of radius a in phase space 
restricted to the finite-horizon regime. 

analytical arguments? In this section we study the idea of Machta & Zwanzig IIMZ83II to approxi- 
mate the deterministic motion by a random walk. 

3.4.1. Derivation of Machta-Zwanzig random walk approximation 

The idea of Machta & Zwanzig IIMZ83II is as follows. If the lattice spacing is such that the discs 
are very close, then the particle will on average be trapped for a long time in each unit cell. We 
refer to such a cell as a trap, where the exits of a trap should have a length which is small compared 
to its total perimeter. Due to the scattering nature of the boundaries, the velocity autocorrelation 
function decays fast, in fact exponentially in the number of collisions IICYOOl . Thus when the 
particle leaves a trap, its motion will be almost uncorrelated with its motion on entering the trap. 
We thus try to approximate the deterministic motion by a completely uncorrelated random walk 
between traps. 

Machta & Zwanzig gave a physical argument which enabled them to calculate the mean res- 
idence time p in a trap. The result agrees with a rigorous calculation discussed below. Once the 
mean residence time is known, we can derive the diffusion coefficient of the random walk via 




(3.26) 
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for a 2D random walk on an isotropic lattice with lattice spacing I and residence time p. 

We remark that the same idea of approximating irregular deterministic dynamics by a stochas- 
tic process has also been explored extensively in the context of transport in area-preserving maps, 
where irregular motion occurs only in part of the phase space: see IIMMP841 Wig92 ILL92I . 



3.4.2. Application to our model 

Let Q be the trap region and ^exitQ the part of its boundary (lying on the edge of the torus viewed 
as a square) which particles can cross, and denote by |A| the m-dimensional Lebesgue measure of 
the m-dimensional set A. Then by an argument detailed in Sec. 13.5.31 together with a method of 
Chernov (see IICM01llChe97l ). we have 



P = — = 73 — ^1 (3-27) 



Here Cy' := (l^exitQI • jS'^ '|) 'is the normalising constant for the measure v' introduced in 
Sec. I3.5.3] and := . |)^' is the normalising constant for the measure jj.. Further, 
B^-i := {x G M"'-i : ||x|| ^ 1} is the (<i- l)-dimensional unitball and5"'"i :={xeM^: ||x|| = 1} 
is its boundary, the {d — 1) -dimensional unit sphere. In 2D, we have l^"^'^' | = 2n and | = 2, 
so that 

P = T$%r- (3.28) 

1 4x112 1 

The exact expression (13.271) for the mean residence time is precisely analogous to the exact expres- 
sion for the mean free path discussed in ||Che97l ; it also agrees with the more physical argument 
of IIMZ83II . 

For the square Lorentz gas that we consider, there are two different regimes in which we can 
expect the Machta-Zwanzig method to be valid: see Fig. 13.131 In (a), the a-discs are almost 
touching, so that the a-a gaps control the escape process from each trap. In (b), the b-disc is so 
large that the a-b gap will now control the dynamics. 

In regime (a) we have 

q(i) =r^-^(a^ + b^)- ^exitQ^'^ =4(r-2a); l = r (3.29) 
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(a) 



(b) 

Figure 3.13: Two regimes where we expect the Machta-Zwanzig approximation to be valid: (a) r close to 

2, so that the a-discs are close to touching; (b) h close to r/\/2 — a, so that the h- and a-discs are close to 
touching. In each case, the shaded region indicates the trap shape used in the Machta-Zwanzig argument. 
In (b), the dashed Unes show the square lattice of traps whose centres are shown by dots. 



and hence 



whilst in regime (b), we have 



(1) = 



r — 2a 



(3.30) 



Q(2) 



= 4 



LV2 



{a + b) 



/ = -=, (3.31) 

\/2 



and hence 



Introducing the quantities 



L>\"^ — 



'ML 



K — ■K{a^ + b'^y 



(3.32) 



wi:=r — 2a; W2 := —j= - {a + b) , 

v2 

which are the the a-a and a-b gap sizes, respectively, we can summarize the above as 

r2 



Wi, i = l,2, 



(3.33) 



(3.34) 



where the ith approximation is expected to be valid when is small relative to the perimeter of a 
cell. 
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3.4.3. Comparison of Machta-Zwanzig approximation with data 

Figure 13.141 and Fig. 13.151 show comparisons of the two Machta-Zwanzig approximations with 
numerical data for the diffusion coefficient as a function of b and of Q, respectively, where d := 
tw^^(b/a) is again the angle of the point in parameter space from the line ^ = 0. We see that the 
Machta-Zwanzig approximations are good when the respective w,- are small, but fail elsewhere. 

To some extent, these approximations explain the observation from Sec. |3.3| that there is a non- 
trivial maximum for values of r near 2, since in that case the first approximation is valid for small 
b and the second is valid for large b, and they combine to predict a maximum at an intermediate 
value of b. It is clear, however, that other effects are important in this intermediate regime. 

For larger values of r, decreasing b has the effect of allowing the particles freer passage through 
the system, so intuition would again predict the observed increase in the diffusion coefficient on 
decreasing b. 

3.4.4. Trap residence time distribution 

Since the above arguments involve the residence time p, we here study numerically its distribution 
as a random variable, where p : M' ^ M. (The space M' and measure v' are discussed in more 
detail in Sec. l3.5l ) To do this, we need to distribute particles uniformly with respect to the invariant 
measure v' on M' , i.e. on the trap exits with inward-pointing velocities. It is sufficient to distribute 
them uniformly with respect to Liouville measure in the billiard domain and evolve them forward 
until the first intersection with the trap exit boundary. By definition of the measure v', they will 
then be correctly distributed. 

We simulate the dynamics until the first exit from this trap and record the trap residence time, 
for an ensemble of N particles. We then estimate the density of this distribution using a histogram. 
The results are shown in Fig. l3.16l for two different geometries. 

An argument similar to that of Machta & Zwanzig, looking at the area available to escape 
in a small time Af, would imply an exponential distribution. However, there is a minimum trap 
residence time given by the time it takes to enter a trap in the centre of one of its sides with 
a velocity perpendicular to that side, collide with a Z^-disc, and exit along the same path. This 
minimum time is thus 2{r/2 — b) =r — 2b. In Fig. lSTTHl we thus compare the numerical distribution 
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Figure 3.15: Comparison of Machta-Zwanzig approximations with numerical diffusion coefficient as a 
function of 9 for different values of the radius a in parameter space. 
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with an exponential distribution with density 



/expW := ;^exp(Vp) 



(3.35) 



with mean 1/A := p, and a delayed exponential distribution given by 



/delayed exp ('"-> Of 



= < 



p-a 



if ^ a 



(3.36) 







ifx < a 



where a:=r = 2b, both of which have mean p. We see that the delayed exponential captures well 
the correct asymptotic behaviour, whereas the standard exponential does not. 

We remark that we can replace the simple Machta-Zwanzig approximation with an a priori 
more general one where we use a continuous -time random walk formulation with a (delayed) 
exponential residence time distribution. However, it turns out that we obtain the same result for 
the diffusion coefficient from any residence time distribution with a finite mean, as was already 
found in IIShl74ll : see also IIWei941 Chap. 3], and Sec. l5.5l for further discussion of continuous-time 
random walk models. 

3.5. Generalised random walk models and Green-Kubo approaches 

Klages & Dellago MKDOOII extended the Machta-Zwanzig approach in a heuristic way, by includ- 
ing the probability to follow specified sequences of traps. Klages & Korabel | KK02 | then found 
a way to include such corrections in an exact Green-Kubo-type expansion. We believe, however, 
that while the result was correct, and the method was mostly correct, there is a gap in the justi- 
fication, since the relationship between the discrete-time and continuous-time approaches studied 
above was not treated correctly. Here we show how to close the gap in the argument. 

3.5.1. Original argument of Klages & Korabel IIKK021 

In the appendix of IIKK02II . Klages & Korabel argued as follows. We wish to generalise the idea 
of a hopping process between traps, for which we must pass from the continuous-time expression 
for the diffusion coefficient D in terms of the Einstein formula for a symmetric system. 




(3.37) 
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Figure 3.16: Trap residence time distribution, (a) is for r — 2.01, b ~ 0.1 and (b) for r = 2.5, b — 0.5. 
(c) and (d) are the same as (a) and (b), respectively, but plotted on semi-log scales. The continuous curves 
show the delayed exponential density given by (13.36b . 
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to a formula in terms of some discrete time n. The mean (•) here is with respect to the invariant 
measure jj. on the phase space Ai of the flow. 

To do this, Klages & Korabel IIKK02II sej] 

x„:=x(?ip), (3.38) 

where x{t) is the random variable giving the position at time t and x„ is the position at the nth 
step. Note that here the '«th step' simply means n times the (so far arbitrarily chosen) time step p. 
Then (1337] ) impUes that 

D=Iim-^(x2), (3.39) 

where the mean is still with respect to /X. 

Since we want an expression in terms of traps, we must show that we can replace x„ by X„, the 
centre of the trap in which x„ lies. Set x„ := x„ — X„, the displacement of the position x„ at step n 
from the centre of the respective trap. Then (13.371 ) implies that 

D=lim-L((X„+Ax„)'), (3.40) 

where Ax„ := x„ — xo> assuming that all initial conditions are in the trap at 0. We expand 

((X„ + Ax„)2) = (X2 + 2X„-Ax„-|-Ax2) (3.41) 

and note that the last term is bounded by r^, where r is the diameter of a unit cell. The Cauchy- 
Schwarz inequahty for the second term then gives 

|(2X„.Ax„)| ^2(x2)V2(Ax2)i/2^,(X2)i/2, (3.42) 
so that the second term grows only as fast as the square root of the first term. Hence we obtain 

D=lim-4(X2), (3.43) 

solely in terms of the trap at time n. 

^In IMZ83I and (KK02 |, t denotes the average residence time in a trap. We follow the convention in the mathemat- 
ical literature (e.g. ICMOl l ). where T is used for the free path function, with mean t. Denoting by p the trap residence 



time function, we replace the z of IKK02I by p. 
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The next step is to obtain a Green-Kubo formula in terms of the jump vector j„ := X„+i — X„ 
using the telescoping sum 

X„ = Xo+£ji = 'j[;j^, (3.44) 

k=0 k=0 

where the equality follows since Xq = 0. This reduces (X^) exactly to the variance of a sum of 
stationary random variables, giving 



{Xl)=Y,(jk-ji)=n(jl)+2Y,in-m){jo-u), (3.45) 

kJ=0 m=\ 



^ = (jo) + I (l - -) (jo -im) = il + I (jo ■im) + - I m (jo J.) . (3.46) 



so that 



m=l " " ^ m=l " m=l 

Thus if (jo • jHi) decays sufficiently fast as m — > oo, then we can take the limit of (13.461 ) as « ^ oo to 
obtain the Green-Kubo relation 

1 1 °° 

^=T^(jo> + ^I(jo-j„,). (3.47) 

Klages & Korabel then use this result to try to extend the Machta-Zwanzig approximation, 
as follows. They follow each trajectory and record the sequence of traps visited via a symbol 
sequence. They numerically calculate the probability of occurrence of each symbol sequence and 
use sequences of length up to n as an ?ith-order approximation to the Green-Kubo expression, 
where the Machta-Zwanzig approximation can be regarded as the Oth-order approximation. They 
obtain results which numerically converge exactly to the (numerically) exact diffusion coefficient 
over the whole finite-horizon regime of the triangular Lorentz gas. The same method was applied 
in a different model in IIHKG02L 

3.5.2. Re-examination of Klages-Korabel derivation 

We argue that the derivation given is not directly related to the numerical method subsequently em- 
ployed, although both are independently correct. The reason for this is confusion in the definition 
of the nth time step: 



(i) In the derivation, the rath time step was taken as looking stroboscopically at the flow at the 
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start of each period of length p. 

(ii) In the numerical method, however, the nth 'time step' instead refers to the time when the 
phase space point crosses into the «th trap. 

Both interpretations correspond to valid reductions from the continuous-time flow to a discrete- 
time map, but they are not equivalent: (Q) corresponds to looking at statistical properties of the 
time-p map of the flow, while we will show that ^ corresponds instead to looking at statistical 
properties of a particular Poincare map, which differs from the usual billiard map and which we 
call the torus-boundary map. 

The part of the derivation after the definition of the position x„ 'at the rath time step', is valid in 
each case, and shows that in any of these cases it is equivalent to study the asymptotic behaviour 
either of x„ or of the corresponding trap centre X„, or of the jump vectors j„ defined as above, 
provided their correlations decay sufficiently fast. As far as we are aware, no rate of decay of 
correlations is yet available for the time-p map (which would be related to the still-unsolved 
problem of the correlation decay rate for the billiard flow), nor for the torus-boundary map (which 
does not seem to have been studied previously). 

3.5.3. Torus-boundary map 

Analogously to the billiard map T (which takes one collision with the scatterers to the next), we 
introduce a map T' , which we call the torus-boundary map, mapping one intersection with the 
trap/torus boundary ^exit Q to the next. Thus the domain of T' is 

M:={(q,v)eA^: qG^exitQ, v-n(q)>0}, (3.48) 

and we follow the trajectory up to the edge of the torus and then continue to the point on the 
opposite edge with which that point is identified. 

The torus-boundary map will be undefined for phase points whose trajectories get trapped 
bouncing between the scatterers, i.e. for points on the stable manifold of the fractal repeller. How- 
ever, this set has measure zero IIGas98l . so that T' is defined almost everywhere. 

We can now view the billiard flow as a suspension flow over the torus-boundary map T' , under 
the trap residence time function p, since specifying a phase space point in M' and a value for the 
residence time function p specifies a unique point of M. . The invariant measure v' on the space 
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M' is defined analogously to the invariant measure v for the billiard map T , i.e. by projecting to 
M' , and then by App. |Bl T' is ergodic. Hence, by the Birkhoff ergodic theorem, for almost every 
X G M', we have that 




(3.49) 



i.e. time averages of the residence time function p are equal to p, so that p is indeed the mean trap 
residence time. 

We can now use the relation between the variance of the suspension flow and the variance of 
the base transformation (see App. E]) to get an expression for D in terms of statistical properties 
of the base transformation. This reduces exactly to ( 13.391 ). where x„ is the position at the «th 
crossing of a trap boundary and the average (•) is with respect to the invariant measure V for the 
torus-boundary map. 

The regularity of the roof function required for this result (see App. |B]| is satisfied, since the 
roof function, which is here given by the trap residence time function, appears to decay expo- 
nentially; this was shown numerically in Sec. 13.4.41 although we are not aware of any rigorous 
results. 

The remainder of the above derivation of Klages-Korabel now goes through to give a Green- 
Kubo formula in terms of the torus-boundary map, provided that correlations decay sufficiently 
fast. As far as we are aware, the only results on decay of correlations for billiards are for the billiard 
map, where exponential decay of correlations is proved IIYou98[ IChe99[ ICYOOl . We expect that 
correlations for the torus-boundary map also decay exponentially fast. However, the proofs for 
the billiard map rely on hyperbohcity due to the scattering nature of billiard boundaries. In the 
case of the torus-boundary map, the boundaries are not scattering: the scattering occurs between, 
rather than at, intersections with the boundary. It is thus possible that a rigorous proof could be 
even harder to obtain than for the billiard map. If we nonetheless assume that these correlations do 
decay sufficiently fast, we do have a Green-Kubo expansion such that the Machta-Zwanzig result 
occurs as the 0th order approximation. This completes the justification of the Klages-Korabel 
method. 

We remark that this method extends to any other map such that we can express the billiard flow 
as a suspension over that map. However, we are not aware of other physically interesting maps for 
which that is the case. 
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3.5.4. Low-order approximations 

If the details of the above justification can be made rigorous, we would then be able to rigorously 
assert that the Machta-Zwanzig approximation is indeed a zeroth-order truncation of the infinite 
Green-Kubo sum. In IIKK02L it was observed that a certain simple low-order expression involving 
colUsionless flights across traps gave a very good approximation of the numerical diffusion coef- 
ficient over the whole finite horizon regime for the triangular Lorentz gas. This was an anomaly 
from the point of view of the expansion used in that paper, since this truncation did not appear 
explicitly in the Green-Kubo sum used there. 

We beUeve that this anomaly can be explained by noting that their expression was instead a 
natural low-order truncation of the discrete-time Green-Kubo formula obtained from the time-p 
map, rather than from the torus-boundary map: there we must include in the first term a sum 
over all accessible traps after a time p , which includes at least part of the probability of collision- 
less flights, together with part of the backscattering probability which also appears in their best 
low-order truncation. Their observation then suggests that the Green-Kubo expression discussed 
above converges more rapidly as the truncation order n increases than that used in their paper. 
This faster rate of convergence depends on correlations for the two respective maps; as such, we 
would not expect to be able to predict which of the two expansions would give better low-order 
approximations without a detailed knowledge of these correlations. 

3.6. Reducing the geometrical symmetry 

In this section we study the effect of altering the geometry of the system to reduce the symmetry 
of the unit cells; again our model is a good candidate for this, whereas the triangular model would 
perhaps be less natural. 

3.6.1. Symmetry of diffusion tensor 

We first consider the effect of symmetry on the diffusion tensor. Suppose we measure the diffusion 
coefficients in the original coordinate system and in a new orthonormal coordinate system x' based 
at the same origin, given by x' = R • x, where R is the orthonormal change of basis matrix (i.e. a 
combination of rotations and reflections). Then the diffusion coefficients with respect to the new 
axes are given by 




^At)) = RikRj,\im^{xk{t)xi{t)) =RikRjiDki = (RDR^),.., 



(3.50) 



3.6. Reducing the geometrical symmetry 



73 



where we have used the Einstein summation convention (sum over repeated indices). If the system 
is symmetric under some symmetry element R, then for that particular R we have D-y = D;y, 
since the system looks the same after the transformation. We remark that this idea is known as 



Neumann's principle in the literature on crystal properties |Nye85|, where it is an empirically 
based fact. In our setting, however, it is a theorem, since Liouville measure is invariant under the 
point group of the lattice. 

Thus symmetry under the transformation R implies that the diffusion tensor satisfies 

RDR"^ = D. (3.51) 

This equation in general restricts the allowed values of the diffusion coefficient. For example, if 
the system is unchanged under reflection in the x-axis, given by the transformation y^—y having 
matrix 

/l o\ 

R = , (3.52) 

then we find that (13.511) implies that Dxy = Dy^ = for a 2D system, so that the diffusion tensor is 
diagonal in this coordinate system. In the case of 3D systems (which we have not investigated), 
tables of crystal symmetries and their consequences for possible forms of the 2nd-order diffusion 
tensor can be found e.g. in 088[[Nye85[. 



3.6.2. Displacing central disc 

Here we consider the effect of displacing the Z7-discs away from the centre of the cell^ , parallel to 
the X-axis. Then the system remains symmetric under reflection in the A-axis, but is now asym- 
metric under reflection in the y-axis. By the above, the diffusion tensor must remain diagonal, but 
we may now have D^x 7^ Ay, which is not possible for the square model. 

Starting from geometrical parameters {a,b) lying in the finite horizon regime, we displace the 
b-disc by a distance g along the x-axis in the direction of increasing x. To continue to block the 
central vertical channel, we cannot move the ^-disc too far: we require ^ + g — b < a, i.e. 

g<a + b-^. (3.53) 



I would like to thank Leonid Bunimovich for suggesting this, with the idea that it could help to understand the 
validity of the Machta-Zwanzig approximation. 
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The Z7-disc touches the a-discs when 



{a + bf 



r \ 2 /r 
2-^) +(2 




(3.54) 



so that the discs are not touching when 




r 



(3.55) 



It is possible that there is a finite horizon even if the b- and a-discs touch. In this case, there is no 
possibility for the particle to escape from the vertical channel in which it begins, so that D^x = 0; 
however, there may still be a non-trivial diffusion coefficient Dyy in the j-direction. Note, however, 
that the analysis of IIBS811 IBSC91I does not immediately apply to this case, since the (convex) 
scatterers are no longer disjoint. 

Figure 13.171 shows numerical data for the diffusion coefficients Dxx and Dyy for this system as 
a function of the displacement g. For ^ = 0, we have Dxx = Dyy by the symmetry results above. 
For non-zero g, however, Dxx decreases and Dyy increases. 

Increasing g past the point where the discs touch, the area available in each cell increases; from 
the figure we see that Dyy suddenly increases its growth rate for two of the three sets of geometrical 
parameters studied. 

On the contrary, for r = 2.01 the value of Dyy seems to stabilise for larger values of g, indicating 
some kind of balance between the increased phase space volume available and the control of the 
dynamics by the small trap exits. 

3.6.3. Rectangular model with less symmetry 

We could also stretch the unit cell to a rectangular one, imposing different lattice spacings in the 
horizontal and vertical directions; this introduces two new geometrical parameters. In particu- 
lar, the standard triangular periodic Lorentz gas can be considered as a particular case of such a 
stretched model. Again the symmetry of the diffusion tensor will be reduced in general, although 
not in the special case of the triangular periodic Lorentz gas. 
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Figure 3.17: The diffusion coefficients D^x and Dyy in the x- and j-directions, respectively, as a function 
of the displacement g of the /j-disc from the centre of the unit cell in the x-direction. (a) r = 2.01, b ~ 0.3. 
(b)r = 2.1,^ = 0.4. (c) r = 2.5, ^7 = 0.6. 



Chapter 4 



Fine structure of distributions and the central limit 
theorem in 2D periodic Lorentz gases 



In the previous chapter we studied the statistical properties of the displacement as a function of 
time in terms of means. Those means are taken over the probability distribution of the observable, 
whose shape we investigate in this chapter. 

As described in Chap. |2l it was proved in IIBS81[ IBSC91II that 2D periodic Lorentz gases 
with finite horizon and disjoint scatterers (with sufficient smoothness of the scatterer boundaries, 
namely piecewise C^) satisfy a central limit theorem: the rescaled displacement distribution con- 
verges in distribution to a limiting normal distribution: 

> z as f ^ oo. (4.1) 

In this chapter we study the structure of position and displacement distributions at finite time t. 
We show that they possess afine structure, consisting of a periodic oscillation superimposed on the 
Gaussian shape that we expect from a diffusive process. We will find an analytical description of 
this fine structure in terms of the geometry of the billiard domain, and provide extensive numerical 
evidence that this is indeed the main influence on the fine structure. This gives a physical picture 
of the weak type of convergence occurring in (14.11) . and leads to a conjecture on a possible stronger 
result; we can also give an intuitive estimate of the rate of convergence to the limiting distribution. 

The results in IIBS811 IBSC91II show how we can smooth away the fine structure to obtain 
rigorous proofs of convergence. Our analysis reinstates the fine structure to give a picture of 
how this convergence occurs, making explicit the obstruction that prevents a stronger form of 
convergence to the limiting normal distribution by showing how density functions fail to converge 
point wise to Gaussian densities. 
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4. Fine structure of distributions and central limit theorem 



4.1. Structure of 2D position and displacement distributions 

4.1.1. Statistical properties of position and displacement distribution 

The displacement distribution occurs naturally in the central limit theorem (Sec. l2.3l ) and in Green- 
Kubo relations l|Dor991 lGas98i . whereas the position distribution is more natural if we are unable 
to track the paths of individual particles. Their statistical properties are very closely related, as 
shown by the following discussion for the ;c-component. 

Expanding the mean squared displacement as 

(Axj) = (xj) - 2{xtXo) + (xl), (4.2) 

and applying the Cauchy-Schwarz inequality to the second term on the right hand side, as in 
IIKK02L which gives 

\{x,xo)\^^{xj){xl), (4.3) 

we see that (Axf) and (xj) have the same asymptotic growth rate, so that they both grow linearly 
if one of them does. 

Further, since |;co| ^ §, we can show that the \/7-rescaled position distribution also satisfies 
the central limit theorem if the displacement distribution does, with the same limiting normal 
distribution. From the point of view of statistical properties it is hence equivalent to study either 
the position or the displacement distribution. 

4.1.2. Shape of 2D distributions 

Figure |4!T] shows scatterplots representing 2D position and displacement distributions for a repre- 
sentative choice of geometrical parameters. Each dot represents one initial condition started in the 
central unit cell and evolved for time t = 50; N = 5 x 10"* samples are shown, started from random 
initial conditions distributed uniformly with respect to Liouville measure in the central unit cell. 
Both distributions show decay away from a maximum in the central cell, an overall circular shape, 
and the occurrence of a periodic fine structure. 

These figures are projections to the unfolded billiard domain Q^nf C M? of the density in the 
unfolded phase space Quni x S^. Since the dynamics on the torus is mixing HCMOIL the phase 
space density converges weakly IILM94I to a uniform density on phase space corresponding to the 
invariant Liouville measure: see App.O Physically, the phase space density on the torus develops 




Figure 4.1: (a) 2D position distribution; (b) 2D displacement distribution, r = 2.5; b = 0.4; t = 50; 
A'^ = 5 X 10"* initial conditions. 
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a complicated layer structure in the stable direction of the dynamics: see e.g. ||Dor99ll . Project- 
ing corresponds to integrating over the velocities; we expect this to eliminate this complicated 
structure and result in some degree of smoothness of the projected densities on the torus, and so 
presumably also of the unfolded projected densities. However, we are not aware of any rigorous 
results in this direction, even for relatively well-understood systems such as the Arnold cat map 



These 2D distributions are difficult to work with, and we instead restrict attention to one- 
dimensional marginal distributions, i.e. projections onto the x-axis, which will also have some 
degree of smoothness. We denote the ID position density at time t and position x G M by ft{x) 
and the displacement density for displacement x by gt{x). We let their respective (cumulative) 
distribution functionj^be Ft{x) and Gt{x), respectively, so that 



and similarly for Gj. (When necessary, we will instead denote displacements by The densities 
show the structure of the distributions more clearly, while the distribution functions are more 
directly related to analytical considerations. 

4.2. Numerical estimation of ID distribution functions and densities 

We wish to estimate numerically the above densities and distribution functions at time t from the 
data points {xt^\... ,Xf^^). The most widely used method in the physics community for estimating 
density functions from numerical data is the histogram IIPTVF92I1 . However, histograms are not 
always appropriate, due to their non-smoothness and dependence on bin width and position of 
bin origin IISil86l . In IIARdV02L for example, the choice of a coarse bin width obscured the fine 
structure of the distributions that we describe in Chap. [51 

We have chosen the following alternative method, which seems to work well in our situation, 
since it is able to deal with strongly peaked densities more easily, although we do not have any 
rigorous results to justify this. We have also checked that histograms and kernel density estimates 
(a generalization of the histogram IISil86l ) give similar results, provided that sufficient care is 
taken with bin widths. 

^Henceforth we use 'distribution function', eliminating the redundant 'cumulative': this seems to be the usual 
terminology in probability theory, e.g. IGS92llFel71l . 



IIDor99L 



X 




(4.4) 



— oo 
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We first calculate the empirical cumulative distribution function ||Sco92[ [Sil86l . defined by 
f /"'""(x) := #{i : x^'^ ^ x} for the position distribution, and analogously for the displacement dis- 
tribution. The estimator F,^"^'^ is the optimal one for the distribution function F, given the data, in 
the sense that there are no other unbiased estimators with smaller variance IISco92[ p. 34]. We find 
that the distribution functions in our models are smooth on a scale larger that that of individual 
data points, where statistical noise dominates. (Here we use 'smooth' in a visual, nontechnical 
sense; this corresponds to some degree of differentiability). We verify that adding more data does 
not qualitatively change this larger-scale structure: with = lO' samples we seem to capture the 
fine structure. 

We now wish to construct the density function /, = dFt / dx. Since the direct numerical deriva- 
tive of fj'^^P is useless due to statistical noise, our procedure is to fit an (interpolating) cubic spline 
to an evenly-spread sample of points from F,''™'', and differentiate the cubic spline to obtain the 
density function at as many points as required. Sampling evenly from Ff^"^'^ automatically uses 
more samples where the data are more highly concentrated, i.e. where the density is larger. 

We must confirm (visually or in a suitable norm) that our spline approximation reproduces 
the fine structure of the distribution function sufficiently well, whilst ignoring the variation due 
to noise on a very small scale. As with any density estimation method, we have thus made an 
assumption of smoothness IISil86ll . The analysis of the fine structure in Sec. 14. II justifies this to 
some extent. 

4.3. Time evolution of ID distributions 

Figure l4!2l shows the time evolution of ID displacement distribution functions and densities for 
certain geometrical parameters, chosen to emphasise the oscillatory structure. Other parameters 
within the finite horizon regime give qualitatively similar behaviour. 

The distribution functions are smooth, but have a step-like structure. Differentiating the spline 
approximations to these distribution functions gives densities which have an underlying Gaussian- 
like shape, modulated by a pronounced fine structure which persists at all times (Fig. l4.2f b)). This 
fine structure is just noticeable in Figs. 4 and 5 of IIARdV02ll . but otherwise does not seem to have 
been reported previously, although in the context of iterated ID maps a fine structure was found, 
the origin of which is related to pruning effects: see e.g. Fig. 3.1 of IIKla96l . We will show that in 
billiards this fine structure can be understood by considering the geometry of the billiard domain. 
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Figure 4.2: (a) Time evolution of displacement distribution functions, (b) Time evolution of displacement 
densities, calculated by numerically differentiating a cubic spline approximation to distribution functions, 
r = 2.1;fe = 0.2. 
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Figure 4.3: Definition of the set H{x). 
4.4. Fine structure of position density 

Since Liouville measure on the torus is invariant, if tiie initial distribution is uniform with respect 
to Liouville measure, then the distribution at any time t is still uniform. Integrating over the 
velocities, the position distribution at time t is hence always uniform with respect to Lebesgue 
measure in the billiard domain Q, which we normalize such that the measure of 2 is 1 . Denote the 
two-dimensional position density on the torus at {x,y) G [0, 1)^ by p^°™^{x,y). Then 



(4.5) 



Here, H{x):={y : {x,y) € Q] is the set of allowed y values for particles with horizontal coordinate 
X (see Fig. 14.31 ). and He is the indicator function of the (one- or two-dimensional) set B, given by 



Mb) 



1, ifbeB 
0, otherwise. 



(4.6) 



Thus for fixed x, p^°"^^{x,y) is independent of y within the available space H{x). 



Now unfold the dynamics onto a 1 -dimensional channel in the ;c-direction, as in Fig. 13.51 and 
consider the torus as the distinguished unit cell at the origin. Fix a vertical line with horizontal 
coordinate x in this cell, and consider its periodic translates x + n along the channel, where ?i E Z. 
Denoting the density there by p^'^^annei (^ + n,);), we have that for all t and for all x and y. 



^channel/ , ,^ 



'\x,y). 



(4.7) 
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We expect that after a sufficiently long time, the distribution within cell n will look like the 
distribution on the torus, modulated by a slowly-varying function of x. In particular, we expect 
that the 2D position density will become asymptotically uniform in y within H{x) at long times. 

We have not been able to prove this, but we have checked by constructing 2D kernel density 
estimates IISil86ll that it seems to be true. A 'sufficiently long' time would be one which is much 
longer than the time scale for the diffusion process to cross one unit cell. It does not, however, 
hold for short times: for example, in Fig. 14.41 we see the development of the 2D density at the 
leading edge. At the earlier time, we see that the 'fluid' has streamed past the small disc without 
reaching the region below it; at the later time, this region is beginning to be occupied. At long 
times, the mixing property of the dynamics will have filled the region approximately uniformly. 

Thus we have approximately 

p^"{x,y) ^ p'°---(x,y)p, W = A (4.8) 

where Pt{x) is the shape of the two-dimensional density distribution as a function of .x G R; we 
expect this to be a slowly-varying function. We use '~' to denote that this relationship holds in 
the long-time limit, for values of x which do not lie in the tails of the distribution. Although this 
breaks down in the tails, the density is in any case small there. 

The ID marginal density that we measure will then be given by 

1 

Mx) = J pf"{x,y)dy ^ Mx)h{x), (4.9) 
>.=o 

where h{x) :=\H{x) \ / the normalized height (Lebesgue measure) of the set H{x) at position 
X (see Fig. 14.31 ). Note that H{x) is not necessarily a connected set. 

Thus the measured density ft {x) is given by the shape p, (x) of the 2D density, modulated by 
fine-scale oscillations due to the geometry of the lattice and described by h{x), which we call the 

fine structure fimction. 

The above argument motivates the (re-)definition of Pt{x) so that = h{x)p,{x) with strict 
equality and for all times, by setting 

P,W:=||^ (4.0) 
We can now view Pt{x) as the density with respect to a new underlying measure hX, where A 
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is 1 -dimensional Lebesgue measure; this measure takes into account the available space, and is 
hence more natural in this problem (see also Sec. IC.4I) . We expect that pt will now describe the 
large-scale shape of the density, at least for long times and x comparatively small. 

The conjecture underlying the above heuristic argument is then that for all a; G M, we have 



or equivalently 



h{x) for ally €H{x), (4.11) 



El 1 for ally GH{x), (4.12) 

where the ratio converges to something which is independent of y. We must use this expression 
since the density at any fixed x tends to as ? ^ oo. In fact, regarding both sides as a function of 
y, we could even expect uniform convergence of the form 

ri ^Crl^'L . (4.13) 

where x is still fixed. 

Figure 14.51 shows the original and demodulated densities /; and p, for a representative choice 
of geometrical parameters. The fine structure in /, is very pronounced, but is eliminated nearly 
completely when demodulated by dividing by the fine structure h, leaving a demodulated density 
pt which is close to the Gaussian density with variance 2Dt (also shown). 

Note that although the density has non-smooth points, which affects the smoothness assump- 
tion in our density estimation procedure described in Sec. 14.21 in practice these points are still 
handled reasonably well. If necessary, we could treat these points more carefully, by suitable 
choices of partition points in that method. 

We estimated the diffusion coefficient D as follows. For r = 2.3 and b = 0.5, using = 10^ 
particles evolved to f = 1000, the best fit line for log(Ax^); against log ? in the region f G [500,1000] 
gives (Ax^) ~ 00003^ which we regard as confirmation of asymptotic linear growth. Following 
HKDOOL we use the slope of log(Ajc^), against t in that region as an estimate of 2D, giving D = 
0.1494 lb 0.0002; the error analysis is as in Chap. [3] (Throughout, we denote by g^i the Gaussian 
density with mean and variance a^, and by N^ji the corresponding normal distribution function.) 
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Figure 4.5: Position density ft exhibiting a pronounced fine structure, together with the demodulated 
slowly-varying function pt — ft{x)/h{x) and a Gaussian with variance 2Dt. The inset shows one period of 
the demodulating fine structure function h. r — 2.3; b = 0.5; t = 50. 
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4.5. Fine structure of displacement density 

We can treat the displacement density similarly, as follows. Let rit{x,y) be the 2D displacement 
density function at time t, so that 

X y 

rit{x,y)dxdy = F{Axt^x,Ayti^y), x,yeR. (4.14) 

(Recall that Axt := Xt — xq.) We define the projected versions t]'^''™"''' and t]'™'"*^ as follows: 

7i^"ix,y):=l^7i,ix,y + n), xeR,ye[0,l), (4.15) 

nel 

^^|"^^{x,y)■.= '£r^^"{x + n,y), x,ye[0,l). (4.16) 

Again we view the torus as the unit cell at the origin where all initial conditions are placed. Note 
that projecting the displacement distribution on M? to the channel or torus gives the same result as 
first projecting and then obtaining the displacement distribution in the reduced geometry. Hence 
the designations as being associated with the channel or torus are appropriate. 

Unlike p^°™^ in the previous section, r]/™"** is not independent of t: for example, for small 
enough t, all displacements increase with time. However, we show that 1^^'°™*^ rapidly approaches 
a distribution which is stationary in time. 

Consider a small ball of initial conditions of positive Liouville measure around a point (x,v). 
Since the system is mixing on the torus, the position distribution at time t corresponding to those 
initial conditions converges as f ^ oo to a distribution which is uniform with respect to Lebesgue 
measure in the bilhard domain Q. The corresponding limiting displacement distribution is hence 
obtained by averaging the displacement of x from all points on the torus. 

Extending this to an initial distribution which is uniform with respect to Liouville measure over 
the whole phase space, we see that the limiting displacement distribution is given by averaging 
displacements of two points in Q, with both points distributed uniformly with respect to Lebesgue 
measure on Q. This limiting distribution we denote by ri^°"^^{x,y), with no t subscript. 

As in the previous section, we expect the 3'-dependence of ri^^'^'^^^^{x + n, •) to be the same, for 
large enough t, as that of ri^°"^^{x, •) for x £ [0, 1). However, t]'°™'*(x, •) is not independent of y, as 
can be seen from Fig. |4.6[ which is a projected version of Fig. 14. H b) to the torus (with different 
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Figure 4.6: Scatterplot of the 2D displacement density rjl°'-'^^{x,y) on the torus for r = 2.1, b = 0.4 and t = 
50. This corresponds to the projection of Fig. l4.ir b) to one unit cell (although the geometrical parameters 
used are different). This density is close to the limiting displacement density t]'™'"\ since the relaxation is 
fast. 



geometrical parameters). We thus expect 

ri^"{x,y)c^ri'°"'^{x,y)U^). (4.17) 



To obtain the ID marginal density gt{x), we integrate with respect to y: 

1 

8t{x)= J rif"{x,y)dyc^^{x)f),{x), (4.18) 

>.=o 

where 

1 

(j){x) := J ri'""''{x,y)dy. (4.19) 

v=0 

Again we now redefine fj so that gt (x) = (x) fjt (x) , with the fine structure of gt {x) being described 
by and the large-scale variation by fl{x), which can be regarded as the density with respect to 
the new measure ^ X taking account of the excluded volume. In the next section we evaluate ^ ix) 
exphcitly. 
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4.5.1. Calculation of .Jc-displacement density 0(x) on torus 

Let {Xi,Yi) and {X2,Y2) be independent random variables, distributed uniformly with respect 
to Lebesgue measure in the billiard domain Q, and let AX := {X2 — X1} G [0,1) be their x- 
displacement, where {•} denotes the fractional part of its argument. Then AX is the sum of two 
independent random variables, so that its density (p is given by the following convolution, which 
correctly takes account of the periodicity of h and <p with period 1 : 

1 

(P{^) = Jh{x)h{x + ^)dx. (4.20) 



This form leads us to expand in Fourier series: 

h{x) = £^(/t)e2'^'*^ = A(0)+2£A(A:)cos27r;tx, (4.21) 

keZ keN 

and similarly for (p , where the Fourier coefficients are defined by 

1 1 
h{k) := J /j(x)e"27ri/:^dx = J /j(x) cos(27rfcc)dx. (4.22) 







The last equality follows from the evenness of h, and shows that h{k) = h{—k), from which the 
second equality in (14.211 ) follows. Fourier transforming (14.201 ) then gives 

^{k)=h{k)h{-k)=h{kf. (4.23) 



Taking the origin in the centre of the disc of radius b (see Fig. l4.3l ). the available space function 
h is given by 

Kx) = ^ {l-2^/b^-x^-2^a^-{\-xf^ (4.24) 

for X G [0, 1/2), and is even and periodic with period 1. (Here we adopt the convention that y/a = 
if a < to avoid writing indicator functions explicitly.) The evaluation of the Fourier coefficients 
of h thus involves integrals of the form 

a 

[ COS zt Va^-t^dt = ^Ji{za), (z / 0) (4.25) 
J 2z 
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where J\ is the first order Bessel function; this equality follows from equation (9.1.20) of HAS 7011 
after a change of variables. 

The Fourier coefficients of h are thus ^(0) = Jq h{x) = \ and, for integer k^O, 



Hl^) = -T7^ {-if aJi{27ia\k\) +bJi{27ib\k\) 



(4.26) 



Note that /q ^{x)dx = 0(0) = h{0)^ = 1, so that is correctly normalized as a density function 
on the torus. 

In Fig. l4.7l we plot partial sums 0m up to m terms of the Fourier series for analogous to (14.211 ). 
We can determine the degree of smoothness of 0, and hence presumably of gt, as follows. The 
asymptotic expansion of Ji (z) for large real z (equation (9.2.1) of IIAS70II ). 

7i(2) ^ cos = ^^^"'^')' ^4.27) 

shows that h{k) = 0{k^^l^) and hence ^{k) = 0{k^^). From the theory of Fourier series (see e.g. 
IIKat04[ Chap. 1]), we hence have that is at least (once continuously differentiable). Thus the 
convolution of h with itself is smoother than the original function, despite the non-differentiable 
points of h. 

We have checked numerically the approach of / r]^°™^ {x,y)Ay to ^ (x) , and it appears to be fast, 
although the rate is difficult to evaluate, since a large number of initial conditions are required for 
the numerically calculated distribution function to approach closely the limiting distribution. 

4.5.2. Structure of displacement distribution 

In Fig. |4.8| we plot the numerically-obtained displacement density gt {x), the fine structure function 
calculated above, and their ratio f]t{x), for a certain choice of geometrical parameters. Again the 
ratio is approximately Gaussian, which confirms that the densities can be regarded as a Gaussian 
shape modulated by the fine structure 0. 

However, if r is close to 2a, then f]t itself develops a type of fine structure: it is nearly constant 
over each unit cell. This is shown in Fig. 14.91 for two different times. We plot both gi and f],, 
rescaled by and compared to a Gaussian of variance 2D. (This scaling is discussed in Sec. 14.61 ) 

This step-like structure of f\t is related to the validity of the Machta-Zwanzig random walk 
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approximation discussed in Sec. 13.41 Having fjt constant across each cell indicates that the distri- 
bution of particles within the billiard domain in each cell is uniform, as is needed for the Machta- 
Zwanzig approximation to work. 

As r increases away from 2a, the exit size of the traps increases, and the Machta-Zwanzig 
argument ceases to give a good approximation. The distribution then ceases to be uniform in 
each cell: see Fig. 14.51 This may be related to the crossover to a Boltzmann regime described in 
llKDOOl . 

4.5.3. Comparison of position and displacement distributions 

For long times, both position and displacement distributions converge to the same limiting normal 
distribution, so we expect the demodulated position and displacement distributions at a large but 
finite time t to be close. This is confirmed in Fig. 14.101 

We can ask if there is in general any advantage in studying one or the other of the position 
and displacement distributions. Since = h{k)^, the Fourier coefficients of decay faster as 
— > oo than those of h. This translates to better regularity (smoothness) of than of h, as we have 
seen in the particular example above, whereo/z G C'/^ but G 

Further, the convolution reduces the amplitude of the oscillations of the fine structure. We can 
see this by taking the first Fourier coefficient alone as an initial indication of the amplitude, putting 

h{x) = 1 +Acos27rx, (4.28) 

although there may be significant contributions to the amplitude from higher Fourier coefficients, 
as can be seen in Fig. 14.71 The first term ensures that /j is a density (per unit length), i.e. that h 
integrates to 1, and A is the amplitude of the oscillation. Then 

1 

0(^)= j h{x)h{x + t,)dx= I + \A^ co^nt,, (4.29) 

so that the amplitude of the oscillations in (the lowest-order Fourier approximation of) is ^A^. If 
A ^ 1 then the channel would have height at some point and no particle could pass. So the ratio 
of the amplitude of to that of /j is ^A < 0.5. This argument leads us to expect that indeed the 



^/ is Holder continuous with exponent Q < a < 1, denoted / e C", if and only if \f{x) — f{y)\ ^K\x — y\°'. 




Figure 4.8: Displacement density gt, with demodulated t], compai-ed to a Gaussian of vaiiance 2D. The 
inset in (a) shows the fine structure function (p for these geometrical parameters, r = 2.1; h = 0.2; t = 50. 
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Figure 4.9: Displacement density gt and demodulated rj,, both rescaled by a/F, at f = 200 and t — 1000, 
compared to a Gaussian of variance 2D. The inset again shows the fine structure function <j). r = 2.01; 
b = 0.l. 
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Figure 4.10: Comparison of the position density /,, the displacement density gt, the demodulated position 
density p? and the demodulated displacement density f)t, for r — 2.1, b = 0.2 and t — 50. 
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amplitude of fine structure oscillations for the displacement density should be significantly less 
than that for the position density. 

The above two results indicate that in numerical investigations it is more useful to concentrate 
on the displacement density, since we expect better performance of density estimation methods 
when the densities are smoother and less oscillatory. 

4.6. Central limit theorem and rate of convergence 

We now discuss the central limit theorem as f — > oo in terms of the fine structure described in the 
previous section. 

4.6.1. Central limit theorem: weak convergence to normal distribution 

The central limit theorem requires us to consider the densities scaled by yjt, so we define the 
rescaled densities 



where the first factor of \/t normalizes the integral of to 1 , giving a probability density. Fig- 
ure |4TI] shows the densities of Fig. |4.2r a) rescaled in this way, compared to a Gaussian density 
with mean and variance 2D. We see that the rescaled densities oscillate within an envelope which 
remains approximately constant, but with an increasing frequency as f ^ oo; they are oscillating 
around the limiting Gaussian, but do not converge to it pointwise. See also Fig. 14.91 

The increasingly rapid oscillations do, however, cancel out when we consider the rescaled 
distribution functions, given by the integral of the rescaled density functions: 



Figure |4~T2] shows the difference between the rescaled distribution functions and the limiting nor- 
mal distribution with mean and variance 2D. We see that the rescaled distribution functions do 
converge to the limiting normal, in fact uniformly, as f ^ oo; we thus have only weak convergence 
of the distributions. 

Although this is the strongest kind of convergence we can obtain for the densities with 
respect to Lebesgue measure, Fig. l4.9l provides evidence for the following conjecture: the rescaled 



gt{x) := VF^f(xV7) 



(4.30) 



X 




(4.31) 
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X 



Figure 4.11: Displacement densities as in Fig. l4~2[ b) after rescaling by \/t, compared to a Gaussian density 
with mean and variance 2D. r — 2.1; b — 0.2. 
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Figure 4.12: Difference between rescaled distribution functions and limiting normal distribution with 
variance 2D. r = 2.1; b = 0.2. 
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densities f}t with respect to the new measure hX converge uniformly to a Gaussian density. This 
characterizes the asymptotic behaviour more precisely than the standard central limit theorem. 

4.6.2. Rate of convergence 

Since the (5, converge uniformly to the limiting normal distribution, we can consider the distance 
\\Gt — A'^2d||^> where we define the uniform norm by 

||F|U:=sup|f(x)|. (4.32) 
We denote by N^ji the normal distribution function with variance a^, given by 

X 

N„,{x):=—^ [ ^^/^^"ds, (4.33) 

S— — 00 

Figure l4~T3] shows a log-log plot of this distance against time, calculated numerically from the 
full distribution functions. We see that the convergence follows a power law 

||G,-A^2D,.|L~?"", (4.34) 

with a fit to the data for r = 2.05 giving a slope a ~ 0.482. The same decay rate is obtained for 
a range of other geometrical parameters, although the quality of the data deteriorates for larger r, 
reflecting the fact that diffusion is faster, so that the distribution spreads further in the same time. 
Since we use the same number A/^ = 10^ of initial conditions, there is a lower resolution near x = 
where, as shown in the next section, the maximum is obtained. 

In llPen021 it was proved rigorously that a ^ ^ ~ 0.167 for any Holder continuous observable 
/. Here we have considered only the particular Holder observable v, but for this function we see 
that the rate of convergence is much faster than the lower bound proved in fPlnOll. 

4.6.3. Analytical estimate of rate of convergence 

We now show how to obtain a simple analytical estimate of the rate of convergence using the fine 
structure calculations in Sec. 14.41 

Since the displacement distribution is symmetric, we have Gt{x = 0) = I /2 for all t. The max- 
imum deviation of Gt from Njo^^t occurs near to ;c = 0, where the density function is furthest from 
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Figure 4.13: Distance of rescaled distribution functions G, from limiting normal distribution N20 in log 
log plot. The straight line is a fit to the large-time decay of the data for r = 2.05. 



a Gaussian, while the fine structure of the density gt means that Gt is increasing there (Fig. I4.12I ). 
Due to the oscillatory nature of the fine structure, this maximum thus occurs at a distance of 1/4 
of the period of oscillation from x = 0. 

Since the displacement density has the form gt{x) = (^{x)f\t{x), after rescahng we have 

U^) = ^{xVt)f\{x), (4.35) 

where f]t{x) := \/tf]t{x\/t) is the rescaled slowly- varying part of gt, and the fine structure at time 
t is given by 

^{x^ft) = l+2Y^^{k)co^,{lKkxVt)■ (4.36) 
The maximum deviation occurs at 1 /4 of the period of (xy^), i.e. at x — that 

l/4v/7 



^ ^{k) cos(27r/cx-v/z')dx 

km 



V'/teN,*:odd ^'^'^ 



(4.37) 
(4.38) 
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The correction due to the curvature of the underlying Gaussian converges to as f ^ oo, since this 
Gaussian is flat at ;c = 0. Hence \\Gt -N\\„ = 0{r'^/'^). 

This calculation shows that the fastest possible convergence is a power law with exponent 
a = 1/2, and provides an intuitive reason why this is the case. This should be compared to the 
rate for convergence of rescaled distribution functions corresponding to solutions of the 
diffusion equation in App. El If the rescaled shape function f\t converges to a Gaussian shape at a 
rate slower than f^'/^, then the overall rate of convergence a could be slower than 1 /2. However, 
the numerical results in Sec. 14.6.21 show that the rate is close to 1/2. We remark that for an 
observable which is not so intimately related to the geometrical structure of the lattice, the fine 
structure will in general be more complicated, and the above argument may no longer hold. 

4.7. Maxwellian velocity distribution 

In this section we consider the effect of a non-constant distribution of particle speedj^. A Maxwellian 
(Gaussian) velocity distribution was used in polygonal and Lorentz channels in IILCW03I and 
IIAACG99II . respectively, in connection with heat conduction studies. The mean squared displace- 
ment was observed to grow asymptotically linearly, but the relationship with the unit speed situa- 
tion was not discussed. 

We show that the mean squared displacement grows asymptotically linearly in time with the 
same diffusion coefficient as for the unit speed case, but that the limiting position distribution may 
be non-Gaussian. For brevity we refer only to the position distribution throughout this section; 
the displacement distribution is similar. 

4.7.1. Mean squared displacement 

Consider a particle located initially at (xo,Vo), where vq has unit speed. Changing the speed of 
the particle does not change the path it follows, but only the distance along the path traveled in a 
given time. Denoting by O^xq, Vq) the billiard flow with speed v starting from xq and with initial 
velocity in the direction of the unit vector vq, we have 

«Dt(xo,vo) = 0'*(xo,vo), (4.39) 



■^I would like to thank Heman Larralde for posing this question, and for the observation that the resulting position 
distribution may no longer be Gaussian. 
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where the flow on the right hand side is the original unit-speed flow. If all speeds are equal to v, 
then the radially symmetric 2D position probability density after a long time t is thus 

giving a radial density 

Pt ir\v) = — ^ exp ( — ^ ) . (4.41) 
' ^ 2Dvt ^ \4Dvt J 

(Throughout this calculation we neglect any fine structure.) 

If we now have a distribution of velocities with density pv{v), then the radial position density 
at distance r is 

oo 

frHr)= J Pr{r\v) pv{v)dv. (4.42) 

v=Q 

The variance of the position distribution is then given by 

oo 

(x2)= J r2/r'(r)dr (4.43) 

r=0 

OO 

= 4Dt j vpv{v)Av=: 4Dtv, (4.44) 


where v is the mean speed, after interchanging the integrals over r and v. 

We thus see that for any speed distribution having a finite mean, the variance of the position 
distribution, and hence the mean squared displacement, grows asymptotically linearly with the 
same diffusion coefficient as for the uniform speed distribution, having normalized such that v = 
1. We have verified this numerically with a Gaussian velocity distribution: the mean squared 
displacement is indistinguishable from the unit speed case even after very short times. 



4.7.2. Gaussian velocity distribution 

Henceforth attention is restricted to the case of a Gaussian velocity distribution. For each initial 
condition, we generate two independent normally-distributed random variables vi and V2 with 
mean and variance 1 using the standard Box-MuUer algorithm IIPTVF92II . and then multiply by 
a, which is a standard deviation calculated below. We use vi and vi as the components of the 
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velocity vector v, whose probability density is hence given by 



p(v) = p(vi , V2) = j= j= = , (4.45) 



where v := |v| = y + is the speed of the particle. The speed v thus has density 

py{v) = ^e-^'y^''' (4.46) 



and mean v = a \/7i/2. To compare with the unit speed distribution we require v = 1 , and hence 
a = y/lpK. As before, we distribute the initial positions uniformly with respect to Lebesgue 
measure in the billiard domain Q. 

4.7.3. Shape of limiting distribution 

The position density (14.421 ) is a function of time. However, the Gaussian assumption used to derive 
that equation is valid in the limit when f ^ oo, so the central limit theorem rescaling 

~f]^\r) := ^tfT\r^t) (4.47) 

eliminates the time dependence in (I4.42I ). giving the following shape for the limiting radial density: 

CO 

~r\r) dv =: ^7, (4.48) 

^ ^ 4D 7 ^ V 4Dv 4 / 4D ' ^ ' 

v=0 



denoting the integral by 7. We can evaluate this integral expUcitly using Mathematica IIWol04ll in 
terms of the Meijer G-functio^ IIEMOT53I IMS73I : 



7 = G, 



3,0/ ?rr4 



(4-51) 

i,0,Oy 



^The Meijer G-function is a generalisation of the classical Gauss hypergeometric function, defined as the following 
Mellin-Barnes-type integral; for more details see IEMOT5 3 1 IMS73] IMKOOl : 



a\,...,ap 
bi,...,bg 



z^'ds; (4.49) 



- m/^ ^7-^- (4.50) 



C is a certain contour in the complex plane and there are restrictions on the a/ and the bj 
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See MMKOOl and references therein for a review of the use of such special functions in anomalous 
diffusion. 

We can, however, obtain an asymptotic approximation to / from its definition as an inte- 

.2 2 

gral, without using any properties of special functions, as follows. Define K{v) := ^ + 
the negative of the argument of the exponential in (14.481 ). Then K has a unique minimum at 
Vmin '■= {r^ / {2TlD)yl^ and we expect the integral to be dominated by the neighborhood of this 
minimum. However, the use of standard asymptotic methods is complicated by the fact that as 
r — > 0, Vmin tends to 0, a boundary of the integration domain. 

To overcome this, we change variables to fix the minimum away from the domain boundaries, 
setting w := v/vmin- Then 

/ = v„,in j ^'-"^(^''^dw, (4.52) 

where a := -^^^^ and L{w) '■= ■^ + ^, with a minimum at Wmin = 1. Laplace's method (see e.g. 
IICKP66II ) can now be applied, giving the asymptotic approximation 

/-Vm.n.-"^(--'-' , ^ , = ^e-'-l\ (4.53) 

valid for large a, i.e. for large r. 
Hence 

/-d(^)'-— Cre-'^'-''', (4.54) 

where 




(4.55) 



4.7.4. Comparison with numerical results 

Figure |4~T4l shows the numerical radial position density ft^'^{r) for a particular choice of geomet- 
rical parameters. We wish to demodulate this as in Sec. 14. 11 to extract the slowly-varying shape 
function, which we can then compare to the analytical calculation. 

The radial fine structure function 0'^''(r) must be calculated numerically, since no analytical 
expression is available. We do this by distributing 10^ points uniformly on a circle of radius r 
and calculating the proportion of points not falling inside any scatterer. This we normalize so that 
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0rad oo, using the fact that when r is large, the density inside the circle of radius r 

converges to the ratio [r^ — n{a^ + b'^)]/r^ of available area per unit cell to total area per unit cell. 
We can then demodulate fj^'^ by 0'^'^, setting 

Figure l4~T4l shows the demodulated radial density Pt^'^{r) at two times compared to the ex- 
act solution (I4.48I) - (I4.51I ). the asymptotic approximation (I4.54I )- (I4.55I ). and the radial Gaussian 
■^e^''^!^^ . The asymptotic approximation agrees well with the exact solution except at the peak, 
while the numerically determined demodulated densities agree with the exact long-time solution 
over the whole range of r. All three differ significantly from the Gaussian, even in the tails. We 
conclude that the radial position distribution is non-Gaussian. A similar calculation could be done 
for the radial displacement distribution, but a numerical integration would be required to evaluate 
the relevant fine structure function. 

An explanation of the non-Gaussian shape comes by considering slow particles which remain 
close to the origin for a long time, and fast particles which can travel further than those with unit 
speed. The combined effect skews the resulting distribution in a way which depends on the relative 
weights of slow and fast particles. 

4.7.5. ID marginal 

The ID marginal in the ;t;-direction is shown in Fig. 14.161 Again there is a significant deviation of 
the demodulated density from a Gaussian. From (I4.54I ). the 2D density at {x,y) is asymptotically 



~£exp[-/3(x2+/)2/3' 



(4.57) 



from which the ID marginal /(x) is obtained by 

oo 

~f{x)= I ~f{x,y)dy. (4.58) 



It does not seem to be possible to perform this integration explicitly for either the asymptotic 
expression (14.571 ) or the corresponding exact solution in terms of the Meijer G-function. Instead 
we perform another asymptotic approximation starting, from the asymptotic expression (14.571 ). 
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3.5 




r 



Figure 4.14: The radial density function f]^'^ compared to the numerically calculated radial fine structure 
function 0'''"^, rescaled to converge to 1 and then vertically shifted for clarity. The demodulated radial 
density p'"'' is also shown, r = 2.3; b = 0.5; f = 100. 




'o 0.5 1 1.5 2 



r 

Figure 4.15: Comparison of the demodulated radial density p,™'' with the exact Meijer-G representation, 
the large-r asymptotic approximation, and the radial Gaussian with variance 2D. 
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Figure 4.16: Rescaled ID marginal of the displacement density g and the demodulated version fj compared 
to the Gaussian with variance 2D and to the asymptotic expression. The latter is not shown close to x = 0, 
where it drops to 0. r = 2.3; b = 0.5. 



Changing variables in (14.581) to z '■= y/x and using the evenness in y gives 

CO 

C\x\ 



m 



Ik 



exp 



(4.59) 



where k := ji \x\'^^^. Laplace's method then gives 



(4.60) 



valid for large x. This is also shown in Fig. 14.161 Due to the \x\ factor, the behaviour near x = 
is wrong, but in the tails there is reasonably good agreement with the numerical results. 



Chapter 5 



Normal and anomalous diffusion in polygonal billiard 

channels 

5.1. Introduction: ergodic and statistical properties in polygonal billiards 

The rigorous results on statistical properties of Lorentz gases discussed in Chap. |2] depend on the 
hyperbolicity of the dynamics due to the defocusing nature of the scatterers. Recently the question 
has been asked to what extent hyperbolicity is actually necessary (rather than just sufficient) for 
these strong statistical properties to hold l|DCv991[DC01ll . 

The weak nature of the scattering from flat obstacles implies that in polygonal billiards the 
Kolmogorov-Sinai entropy and all Lyapunov exponents are 0, so that infinitesimally separated 
trajectories separate slower than exponentially; in fact they separate linearly, as can be seen by un- 
folding the polygonal billard and considering two initial conditions with slightly different angles. 
Nonetheless, trajectories separated by a finite distance may fall on different sides of a polygonal 
vertex, causing them to separate faster. Characterising this effect is difficult: see e.g. llvB04ll for a 
recent attempt. 

If all angles are rational, i.e. rational multiples of n, then the dynamics reduces to directional 
flows on invariant surfaces and is well-understood; we do not consider this case. Rigorous results 
on ergodic properties of more general polygonal billiards are reviewed in IIGut861 lGut96[ lGut03l . 
There are known rigorous examples of ergodic polygonal billiards, but no examples known to be 
mixing, although mixing has not been dispioved in general. Recently numerical evidence has been 
given that the billiard dynamics in right-angled triangles with irrational acute angles is ergodic and 
weak-mixing but not mixing IIACG97i . while in triangles with all angles irrational the dynamics 
appears to be mixing IICP99I . 

Despite these uncertainties, various polygonal billiard models have been found numerically to 
exhibit normal diffusion, in the sense that the mean squared displacement grows asymptotically 
linearly in time (property (a) from Chap.© HDCOI [ lARdV02l ILCW03]| . In fliis chapter we furflier 
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explore statistical properties of polygonal billiards, attempting to understand when we can hope to 
see normal diffusion and to what extent stronger properties, such as the central limit theorem, are 
satisfied. We also investigate when normal diffusion fails and characterise the anomalous diffusion 
that results. 

Recently there has also been much interest in finding simple model systems which show nor- 
mal heat conduction, in the sense that Fourier's law holds. Several simple periodic billiard mod- 
els were found which show this behaviour when placed between thermal reservoirs at different 
temperatures IIAACG99[ IARdV02l ILCW03I . Evidence has been given that the heat conduction 
properties of ID can be predicted once the diffusive properties are known IILW031 IDKU03I . so 
that it is important to characterise the diffusive behaviour. 

We remark that the weak instability in polygonal billiards means that numerical simulations 
should be more reliable than in scattering billiards, in the sense that computed trajectories should 
lie close to true trajectories for a long time; see e.g. MAGROOII . For this reason we believe that 
numerical experiments on statistical properties of polygonal billiards can provide useful informa- 
tion. 

We begin by reviewing existing models which numerically exhibit asymptotic linear growth 
of the mean squared displacement. We then attempt to identify geometrical features of polygonal 
billiards which allow or prevent the occurrence of normal diffusion, and construct two more classes 
of models which have normal diffusion except when a particular geometrical condition is satisfied, 
namely that parallel scatterers exist, when we find anomalous diffusion. 

5.2. Polygonal models exhibiting normal diffusion 
5.2.1. Previous models 

Model of Alonso et al. Alonso et al. IIARdV02ll studied the geometry shown in Fig. |5.1(a)| and 



Fig. |5.1(b)| We fix the angles 0i and <^ and choose d such that the bottom triangles are half the 
width of the top triangle. This determines the ratio of hi to /12 in terms of the angles 0i and ^2- 
We then require the inward-pointing vertices of each triangle to lie on the same horizontal line to 
prevent infinite horizon trajectories. Taking hi+h2 = h = \ then gives 
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(b) 

Figure 5.1: (a) Geometry of the polygonal billiard unit cell of l|ARdV02| . shown to scale with 02 = ^/(2e). 
(b) Part of the polygonal channel with the same parameters. 

with h\ = fiftan^i and = (<i/2)tan(/>2. We remark that in IIARdV02ll it was stated that the area 
\Q\ =dh of the billiard domain is independent of when is fixed. This is not correct, however, 
since by (I5.1I ). d depends on ^2- 

In IIARdV02ll . the parameters <pi = 7i{V5 - l)/8 and 02 = n/q, q eN,3 <q <9 were used. 
For g ^ 5, the mean squared displacement was found to grow like with a in the range 1 to 1 .08, 
indicating normal diffusion. 

For q = 3,4, however, anomalous diffusion was found, for which {Ax^)t ~ t" with a 7^ 1. As 
far as we are aware, there is as yet no explanation for this observed anomalous behaviour, although 
presumably number-theoretic properties of the angles are relevant; in a second paper ||ARdV04| . 
Alonso et al. state that in these cases they find large classes of periodic orbits which are either 
trapped in one cell in the case q = A, resulting in sub-diffusion, or are propagating orbits, resulting 
in super-diffusion for ^ = 3. 

These results appear to be related to the rational angles used. We have instead used a value 
of 02 which is not rationally related to 0i, namely 02 = 7r/(2e) ~ 7r/5.44, where e is the base of 
natural logarithms. In this case we find that (Ax^)f ~ which we regard as asymptotically 

hnear, with D = 0.3796 ± 0.0009. 

Model of Li et al. Figure l5!2l shows the unit cell of the model introduced by Li et al. IILCW03I . 
Unlike the other models we study, this one generates a two-dimensional polygonal billiard with 
finite horizon when unfolded: presumably this was the reason for introducing this more subtle 
shape. 
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2d 




Figure 5.2: Unit cell of the model introduced in PLCWOSl, with 2d = 3, /i = 1.8, w 2.19, = 1 and 
9 ~ {V2— l)7t/2. Note that ^ — lis a transcendental multiple of 7t, whereas 9 is only a quadratic irrational 
multiple of n. 



5.2.2. Necessary conditions for normal diffusion 

Based on the features of the models presented above and those introduced below, and on ex- 
periments with further models, we arrive at the following heuristic ingredients for constructing 
polygonal billiards with normal diffusion: 

(i) avoid vertex angles which are rationally related to n; 

(ii) avoid infinite horizon trajectories; and 

(iii) avoid parallel scatterers. 

As discussed above, point (i) has to some extent a rigorous justification. Point (ii) is related 
to heuristic arguments discussed in Chap. [6] which show that infinite horizon trajectories lead to 
at least weak anomalous diffusion, where the mean squared displacement grows like tlogt. Point 
(iii) is the main observation that we wish to emphasise. As we discuss in more detail in Sec. 15.41 
parallel scatterers seem to cause a channelling effect, resulting in long laminar stretches and thus 
anomalous diffusion. 



5.2.3. Constructing new models with normal diffusion 

Following the above heuristic rules, we now construct two classes of models which seem to exhibit 
normal diffusion for most parameter values, but anomalous diffusion for certain special geometri- 
cal configurations. These new models provide clear evidence of point (iii) in particular. We follow 
IIARdV02l in using strictly ID channels, i.e. ones which cannot be unfolded in the y direction, 
which enables us more easily to ensure a finite horizon. 
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d 




(b) 



Figure 5.3: (a) One unit cell, and (b) several cells forming a ID channel, for a polygonal model derived 
from the Lorentz gas considered in Chap. [3] There is a crossover from normal to anomalous diffusion when 
hi = h. 

Polygonal Lorentz model Our first model is a polygonalised version of the model introduced 
in Chap. [3l restricted to be strictly one-dimensional: we create a polygonal channel and then add 
an extra scatterer at the centre of each unit cell to ensure a finite horizon by blocking the central 
corridor: see Fig. 15.31 

We have defined this model in terms of heights of the scatterers. The angles between the 
scatterers are then given as arctangents of ratios of these heights. In general these angles will be 
irrational multiples of n, although since in our calculations we use simple heights, we will always 
end up with angles which are nongeneric in the sense that they are arctangents of rationals. 

To generate initial conditions uniformly distributed with respect to Lebesgue measure in the 
available space inside the polygon, we use the 'rejection method' IIPTVF921I . as follows. We 
generate points in a rectangle which contains the polygonal billiard domain and reject points which 
do not lie inside it, using the following algorithm ||0'R98I : a point lies inside the polygon if a ray 
extending from the point to the exterior of the polygon crosses the boundary of the polygon an 
odd number of times. 

Simplified version: zigzag model We can simplify the above model by eliminating the central 
additional scatterer. In order still to be able to block horizontal trajectories in the centre of the 
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(a) 




(b) 

Figure 5.4: (a) One unit cell, and (b) several unit cells forming a ID channel, for a zigzag model which 
is a modification of Fig. 15.31 This model also shows a crossover from normal to anomalous diffusion when 
hi = hj. 

channel, we flip the bottom line of scatterers to point up instead of down, resulting in the modefl 
shown in Fig. [ 



5.3. Normal diffusion 

We present evidence that the models discussed above do indeed exhibit normal diffusion. It is 
difficult to distinguish between asymptotic t and tlogt behaviour by examining only the mean 
squared displacemeno so that we look at several different indicators. (Of course, numerical ev- 
idence alone can never be conclusive without a rigorous basis, especially as regards asympotic 
behaviour; nonetheless, the available numerical evidence presented below leads to this conclu- 
sion.) 



5.3.1. Moments 

Following IIAHO03II . we consider the moments Several of these moments are shown as 

a function of time in Fig. 15.51 for the polygonal Lorentz channel with values of the geometrical 
parameters for which we expect normal diffusion, i.e. without parallel scatterers. 

Since < t^, we expect these moments to have power law growth, possibly with e.g. loga- 
rithmic corrections. We thus define the growth rates jq by IIAHO03I1 



lim 



log? 



(5.2) 



^[The same model was studied independently in IIJR06I at around the same time, with similar conclusions. 
■^[This was done in the published version of this chapter 1SL06I .1 
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Figure 5.5: Moments as a function of f for different values of q. Polygonal Lorentz channel with 

hi = 0.1, ho, = 0.45, A/i = 0.1, w = 0.2. 



Corrections to power law growth are ignored by this definition. Numerically we calculate jq by 
fitting a straight line to a log-log plot of the ^h moment in the long-time regime. The growth 
rates are shown in Fig. l5.6l for the polygonal Lorentz channel in the normal diffusion case. We see 
a qualitative change at around q = 3 between two approximately linear regimes. Results of this 
type have been reported in a variety of models in MCMMG V99[ IFM YO 1 [ I AC041 and are sometimes 
referred to as 'phase transitions'. Physically this corresponds to a change in importance between 
ballistic trajectories and diffusive trajectories MFMYOlll . 

We also remark that it was proved in IIAHO03II that jq is a convex function of q, assuming 
that the limit defining jq exists. This considerably restricts the possible behaviour. Again the data 
points on the graph violate this and so cannot be correct; nonetheless we believe that the data 
gives an impression of the true behaviour. In particular, the numerical accuracy of the high-order 
moments is not good, since they become very large: the rate of increase in the high-order regime 
cannot in fact exceed 1, contradicting the best-fit linj^. 



[Higher moments are dominated by extreme values, which must be ehminated to obtain reproducible results 

ISLoel .i 
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Figure 5.6: Growth rate Yg of the ^th moment {\x\'') for the polygonal Lorentz channel with hi = 0.1, 
/z3 — 0.45, A/j = 0. 1, w = 0.2. The straight lines are best fits to the low- and high-order moments. 

5.3.2. Velocity autocorrelation function 

The above results provide some evidence of normal diffusion: we have (x^), so that the 

growth of the variance of the position distribution (and hence of the mean squared displacement) 
is asymptotically linear. We can attempt to confirm this by studying the velocity autocorrelation 
function (vq v^), which must be integrable for the diffusion coefficient to exist. As seen in Fig. 15.71 
this is an oscillatory function whose maxima seem to decay approximately as (Statistical 
errors (not shown) dominate for large time.) Since there is also cancellation due to the oscillation, 
this gives weak evidence that the velocity autocorrelation function is integrable. 

It was suggested in IILM93I to consider instead the integrated velocity autocorrelation function 



since the delicate cancellations in C{t) will be seen more robustly in R(t). If C{t) is integrable then 
R{t) 2D as f ^ oo. If C{t) decays as, for example, C{t) ~ t^^ then R{t) diverges logarithmically, 
so that plotting R{t) against log? should give a linear asymptotic growth. Figure [5^ provides some 
evidence that C{t) is integrable, since an asymptotic limit seems to be attained as t increases. 




(5.3) 
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Figure 5.7: Log-log plot of velocity autocorrelation function (vqv,) as a function of time. The function 
oscillates about 0, so different colours are used to depict the positive and negative parts. A straight line with 
slope —1.05 is shown. 



Thus overall we have good evidence that the mean squared displacement can grow asymptotically 
linearly, although we have not ruled out a small logarithmic correctiorc . 



5.3.3. Fine structure 

The shape of the displacement density was considered in IIARdV02l using histograms, but the 
results were not conclusive. Here we use the more refined methods developed in Chap. |4] to 
describe the fine structure and to show that the position and displacement distributions do seem to 
be asymptotically normal. 

Figure [S!9l shows the position density ft (x) in the Alonso model at a particular time. Following 
the method of Sec. l4.4[ we can calculate the fine structure function h{x) as the normalized height 
of available space at position x. Taking the origin in the centre of the unit cell in Fig. |5.1(a)[ we 
have 



h{x) 



2d 

\Q\ 



[x tan 01 + |x — 2<3? I tan ^2) 



(5.4) 



for ^ ^ (i, with h being an even function and having period 2d, and \ Q\ = dh. (The factor of 
2d makes h a density per unit length.) This fine structure function is shown in the inset of Fig. 15.91 



*[This point is again addressed in ISL06I .1 
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Figure 5.8: Integral (vqAx,) of the velocity autocorrelation function as a function of logj^f. 



We demodulate ft by dividing by h to give Pt{x) := ft{x)/h{x), which is also shown in the figure. 
We see that it is close to the Gaussian with variance 2Dt. 



With the same notation as in Sec. I4.5.1I we can also calculate the fine structure function ^ of 
the displacement density for the polygonal channel. The Fourier coefficients are /j(0) = 1 and 



^W = ^/K-)cos(^ 



1 



m 



for k^O, where for m G Z we have 



4tan((^i), if/: odd 
8tan(02), if A: = 4m + 2 
0, 'ifk = 4m. 



(5.5) 



(5.6) 



For the polygonal Lorentz channel, we have for x G [0,(i] that 



Id 

\Q\ 



2hi+h2 + h- ^^^J^^ x-{w-x)t\f:^„\{x) 



(5.7) 
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where 



This gives 



\Q\=2[2hid + ^dih2 + h)-w^]. (5.8) 



1 — cos Tik) — I 1 — cos 



(5.9) 



d \ d 

In both cases we have h{k) = 0{k^^) and hence ^{k) = 0{k^^), so that is at least C^, whereas 
h is Lipschitz continuous (i.e. Holder with exponent a = 1). 

5.3.4. Central limit theorem 

As for the Lorentz gas, we rescale the densities and distribution functions by to study the 
convergence to a possible limiting distribution. Again we find oscillation on a finer and finer scale 
and weak convergence to a normal distribution: see Fig. l5.10l Figure lSTTTl shows the time evolution 
of the demodulated densities f]t. There is an unexpected peak in the densities near x = for small 
times, indicating some kind of trapping effect; this appears to relax in the long time limit to a 
Gaussian density. Again we conjecture that we have uniform convergence of these demodulated 
densities to a Gaussian density. 

Figure I5TT2] shows the distance of the rescaled distribution functions from the limiting normal 
distribution, analogously to Fig. I4.13[ for several values of 02 for which the mean squared dis- 
placement is asymptotically linear. The straight line fitted to the graph for <^2 = {2e) has slope 
—0.212, so that the rate of convergence for this polygonal model is substantially slower than that 
for the Lorentz gas, presumably due to the much slower rate of mixing in this system. A similar 
rate of decay is found for 02 = ?r/7, while 02 = 6 and 02 = 9 appear to have a slower decay rate. 
Nonetheless, the distance does appear to converge to for all these values of 02, providing evi- 
dence that the distributions are indeed asymptotically normal, i.e. that the central limit theorem is 
indeed satisfied. A similar calculation for the polygonal Lorentz channel (for the parameters used 
above) gives a convergence rate of t^^-^^ , of the same magnitude as for the Alonso model. 

5.4. Anomalous diffusion 

In dispersing Lorentz gases, there is a weak form of anomalous diffusion, with the mean squared 
displacement growing as (Ax^), ~ tXogt for large t: see Chap. [6] In polygonal billiards of many 
types, however, a strong type of anomalous diffusion has been observed, with (Ax^); ^ for some 
a > 1. Note that the maximum possible rate is a = 2, since \x\ < t due to the finite particle speed. 
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Figure 5.10: Rescaled displacement denstities compared to the Gaussian with variance 2D. The inset 
shows the function ^ for this geometry. 
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Figure 5.12: Distance of the rescaled distribution functions from the hmiting normal distribution for the 
polygonal model with different values of 02- The straight line is a fit to the large-time decay of the irrational 
case ^^n/ (2e). 
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Figure 5.13: mean squared displacement as a function of time t, on a log-log plot. h\ — 0.1, h2=hj,= 0.3, 
w — 0.2, = 1, for the polygonal Lorentz channel. The straight line is a fit to the long-time slope. 

As stated above, we find anomalous diffusion in our models precisely when there are parallel 
scatterers in the unit cell. Previous observations of anomalous diffusion in billiard models include 
||ZE97[ IZas02[ ILWH02I . and MPZOU in discrete time, but they do not seem to have explicitly 
related the occurrence of anomalous diffusion to the geometry of the system. We remark that in 
IIZwa83i a model was presented for which the unit cell is a rectangle with a small window in one 
edge. It was stated that the growth exponent of the mean squared displacement depends on the 
number-theoretic properties of the aspect ratio of the rectangle. It would be interesting to revisit 
this model in the light of our results. 

5.4.1. Moments 

Figure [5TT3] shows a typical plot of the mean squared displacement when h2 = h^. The long-time 
behaviour is well described by 

(Ax2),~Cf«, (5.10) 

where a = 1 .40. 

Figure [57T4] shows the growth exponent of the higher order moments {\x\'^)t for the polygonal 
Lorentz channel. There are two different linear regimes, as previously, with a crossover occurring 
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Figure 5.14: Growth exponent of higher order moments for the polygonal Lorentz channel with anomalous 
diffusion. There is again a crossover between two different linear regimes occurring close to ^ = 3. 

near q = 3. Figure |5T5] shows the same for the zigzag model, for which there does not seem to be 
a crossover, and we have strong anomalous diffusion under the definition of IICMMGV99L where 
Yq = vq for all q, with the constant V given by V ~ 0.93 for these parameters. 

5.4.2. ID densities 

Since the mean squared displacement grows like {Ax^)t ~ t" with a > 1, the \/7 scaling used in 
the central limit theorem cannot give densities which converge in any sense, since the variance of 
the rescaled densities would still tend to infinity. Instead we rescale to keep the variance bounded, 
setting 

f,(x):=?«/2g,(x^«/2). (5.11) 

This scaling was previously used in MPZOll for the Poincare map of an anomalously diffusive 
billiard. 

Figure 15.161 shows the .x-position density. The fine structure is again mostly removed by de- 
modulating by the fine structure function h found above, although this demodulation does not seem 
to be as effective as previously; this could be an artifact of an insufficiently precise estimation of 
the original density, or it could reflect weaker mixing properties of the system. 
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fS- 3 




Figure 5.15: Growth exponents 7^ of higher order moments for zigzag channel. The growth exponent is 
close to a linear function of q, although perhaps two different linear regimes can be detected for low and 
high values of q. 




Figure 5.16: Position density and demodulated position density in the polygonal Lorentz model with 
hi =0.1,/z2 = /!3 =0.45, w = 0.2,ci = 1 andr = 100. 
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Figure ISTTTl shows a sequence of demodulated densities rescaled as above. They appear to con- 
verge at long times to a limiting shape which is non-Gaussian. They are compared to a Gaussian 
with variance B, where B is the generalisation of the diffusion coefficient given by (Ax^); ~ Bt". 
The demodulated densities are, however, rather noisy; those for the displacement distribution are 
even more so. 

Zigzag model Rescaled densities for the zigzag model are shown in Fig. 15.181 for long time 
we again have convergence to a limiting distribution. The spikes in the tails seem to correspond 
to long-lived propagating orbits. Note that there is no fine structure to take into account for the 
zigzag channel, since the channel width is a constant within the unit cell. The peak at ;c = has 
been observed previously in anomalous diffusion: see e.g. MFMYOlll . 

We remark that the unit cell of the zigzag model when the top and bottom scatterers are parallel 
can be reduced to a parallelogram with irrational angles. The mixing properties of the billiard flow 
inside such a geometry are not obvious, although we have checked that velocity autocorrelations 
decay as a power law. There may be an effect of slow ergodicity, where the system takes a long 
time to explore certain regions of phase space: see e.g. IIKH98[|PZ01II . 

5.4.3. Maxwellian velocity distribution 

Just as for the Lorentz gas, we can consider the effect of a different distribution of velocities; this 
has been considered for heat conduction BLC W03 [ I ARdV02l . 

Recall from Sec. l4.7l that if there is a speed distribution py (v), then the mean squared displace- 
ment at time t in 2D is given by 



If a billiard system shows anomalous diffusion for particles with unit speed vq = 1 which looks 
hke 




(5.12) 



V 



{r\t)) 



a 



(5.13) 



then 



(r2(0)v = (r^av/vo))vo = (fv) 



(5.14) 
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Figure 5.17: (a) Rescaled demodulated densities in the polygonal Lorentz model for several times, com- 
pared to a Gaussian with the same variance; the central region is shown in more detail in (b). There seems 
to be convergence at longer times to a limiting shape which is non-Gaussian. 
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Figure 5.18: Rescaled densities for the zigzag channel, where (Ax^)f ^ f '-^^ 

Thus 

(Ax^) = ?«y";jy(v)v"dv. (5.15) 

V 

Hence there is now a correction factor of the ath moment of the speed distribution which enters, 
although the rate of growth is still f". We have observed such a correction numerically. 

5.5. Continuous-time random walk model for anomalous diffusion 

A widely-used framework for understanding anomalous diffusion processes is the theory of continuous- 
time random walks (CTRW): see e.g. ||Wei941lKBS871IZK93i This is a generalisation of standard 
discrete-time random walks where individual steps are still independent, but are now described by 
a density function i//^(r,f), where i//^(r,f)drdf is the probability of having a step with distance in 
(r,r+ dr) which takes time m{t,t + At). 

To be able to model super-diffusion, \\f must have a coupled form IIZK93II . such as 

W{v,t) = \5{\r\-t)\if{t). (5.16) 

This form (called the velocity model in IIZK931 ) models motion at a constant velocity for a time t 
in the direction r; after each stretch the direction is randomised and a new step is taken. It gives 
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the following long-time growth of the mean squared displacement o-^{t): 



0<v^l 



t 



3-v 



l<v<2 




(5.17) 



flnf 



v = 2 



t 



V > 2; 



when \^{t) ~ t 



i-v. 



; see IIGei95[|KZS95H for reviews. More generally we have 



\\r(r,t)=p(x\t)\\f(t) 



(5.18) 



where p(^\t) is the conditional probability of a step of length r if we know that it lasts for time 
t: see e.g. ||Wei941 Sec. 2.5] and ISK89II . In general the rate of growth of the mean squared 
displacement will then also depend on the properties of p(x\t). 

Figure 15.191 shows a representative trajectory in the zigzag model with anomalous diffusion. 
From the figure we can see the importance of laminar motion, i.e. coherent motion in one direction 
along the channel. Referring to the above CTRW picture, we try to model the motion using random 
walks where the steps are stretches of laminar motion. 

We define laminar motion in the zigzag model as follows. For each section (half a unit cell) 
/ of the channel, we assign a vector L, parallel to the channel walls, with a consistent ;c-direction 
along the channel. We say that two consecutive stretches of trajectory (between bounces) are part 
of the same laminar motion if they have the same directions along the channel, as measured by 
sgn(v • L,). Figure [5^20l shows a scatteiplot representing the joint distribution \\f{x^t) of periods of 
laminar motion with A;-displacement x along the channel and taking time t. Part (b) of that figure 
shows that allowed values of x are restricted to being near integers, since a laminar motion cannot 
end in the middle of a unit cell. (A similar effect in infinite horizon 2D periodic Lorentz gases was 
described in ||Ble92|| . although no figure was given.) We can thus regard the walk as taking place 
on a one-dimensional lattice. 

We see that the distribution ^f(x^t) is non-trivial, and the figure does not give much hope of 
finding an analytical expression for it. However we could hope that an approximation of the 
form V'(r,f) = ^5(|r| — vt) Yi^) referred to above may be adequate, for some speed v, since the 
distribution is somewhat concentrated along diagonal hnes. Figure \52T\ shows the tail region of 
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Figure 5.20: Scatterplot representing the joint density yf{r,t) of laminar paths in the zigzag model with 
anomalous diffusion {hi = 0.1, /12 = h^ = 0.3). The straight lines have slope 1, representing the maximum 
possible speed, (b) Shows the fine structure near the origin: the allowed values of x are restricted to be near 
integers. 
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Figure 5.21: (a) Density y/{t) of times of laminar motion for the same model as in Fig. 15.201 

the density and its distribution function. If we assume tlie velocity model form, then the 
observed long-time power law decay Y{t) ~ t^^-^^, shown in Fig. 15.211 gives v = 1.58 and hence 
a mean squared displacement <J^{t) ~ f3-i.58 _ ^1.42 jj^j^ compares with the observed growth 
G^{t) t^-^^. The agreement leaves room for improvement, showing that we must use a more 
general CTRW model incorporating information on the complete or perhaps reject any 

CTRW model which omits important information on correlations between consecutive laminar 
trajectories; such correlations can be seen in Fig. |5.19[ for example. 

Figure [5^22] shows a trajectory in the polygonal Lorentz gas model. Despite the fact that anoma- 
lous diffusion is also found in this model when there are parallel scatterers, the mechanism is much 
less clear, since it does not seem possible to identify an obvious candidate 'laminar motion' in this 
model by looking at sample trajectories. 

5.6. Crossover from normal to anomalous diffusion 

Since anomalous diffusion occurs under certain geometrical conditions, whereas normal diffusion 
seems to be more general, it is of interest to ask how the crossover from one to the other occurs 
when the geometry is changed to approach one where anomalous diffusion occurs. We study the 
two models introduced above to show numerically how the crossover occurs. 
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Figure 5.22: Part of a single trajectory of the polygonal Lorentz model with h\ = 0.1, h2 = = 0.3. 
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Figure 5.23: Mean squared displacement as a function of time for different values of Ah = h2 — tending 
to 0, for the polygonal Lorentz model with h\ ^QA, = 0.45 and w = 0.2. Values of Ah shown are, from 
top to bottom, A/i = 0, 10"", 10"^", 10"'. 

In the following we fix all geometrical parameters except for h2, which we vary according to 
/i2 = /J3 + Ah. Anomalous diffusion is found for Ah = 0, with asymptotically normal diffusion 
found for any Ah ^ 0. It is thus of interest to study how the statistical behaviour changes from one 
asymptotic regime to the other as the symmetrical configuration is approached when Ah 0. 

Figure 15.231 shows the mean squared displacement as a function of time for several values of 
Ah. The same data is plotted on log-log axes in Fig. I5.24f a) and as deviations from a linear fit to 
one of the curves in Fig. l5.24f b). We see that as Ah —>■ 0, the curves follow the anomalous diffusion 
curve (with slope > 1) for longer times, before a crossover occurs to asymptotic linear behaviour. 

The intercept of the asymptotic linear growth for A/j / is related to the diffusion coefficient 
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(a) 




(b) 



Figure 5.24: (a) Log-log plot of the mean squared displacement as a function of time, (b) Deviations of 
the log-log plot from a straight line fitted to the long-time part of the lowest curve. As Mi approaches 0, the 
mean squared displacement curve follows that for the anomalous case A/i — for longer and longer times. 
Values of Ah shown are as in the previous figure. 
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-logio(|A/i|) 

Figure 5.25: Log-log plot of the diffusion coefficient D as a function of A/i for the polygonal Lorentz 
model. The hnear growth indicates power law behaviour. 

D{Ah). As Ah — > 0, the intercept moves up, corresponding to an increase in D. Figure ls!25] shows 
the diffusion coefficient in this asymptotic linear regime as a function of Ah, obtained via the 
slope of the mean squared displacement. Note that it is clear from Fig. 15.241 that the asymptotic 
linear growth regime has not yet been reached for the smallest values of Ah, so that the diffusion 
coefficients for those values are expected to be underestimates. 

We obtain a straight line on a log-log plot, and a fit in the region 3 ^ — logjQ(A/i) gives 

logioZ)(Aft) ~ -1.10-0.1521ogio(A/i)). (5.19) 

The diffusion coefficients obtained for negative values of Ah are also shown for h = 0.45, and it is 
clear that the growth rate is the same. We hence obtain the power law behaviour 

D{Ah) ~ 0.08 |A/jr°-^^ as \Ah\ 0. (5.20) 



Note that this agrees with D = oo when Ah = 0. 

We could also consider the crossover time, i.e. the time required to switch from the anomalous 
diffusion regime to the linear growth regime. One possible definition of this time could be as the 
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intersection of two straight line fits: a fit to tiie initial anomalous growth and a fit to the asymptotic 
normal growth. A construction of this type was used e.g. in IICGvB04ll in a different context. 
From Fig. 15.241 we see that this crossover time tends to oo as A/i ^ 0, and we could studjf] the 
dependence of this time on A/i. 

5.6.1. Zigzag model 

There is a similar crossover to anomalous diffusion when h2 = in the zigzag model; the growth 
of the moments is very similar to the polygonal Lorentz channel. Figure |S!26] shows the growth of 
the diffusion coefficients as A/j — > 0. A fit gives 

logio£>(A/j) = -0.320 logio(A/j)- 0.191, (5.21) 

and hence 

D{Ah) ~ 0.644 (A/j)"°-^2 as Ah 0, (5.22) 

so that we again obtain a power law, but with a different rate of growth than in the polygonal 
Lorentz model. 

5.6.2. Qualitative explanation 

In the zigzag model we can find a qualitative explanation of the above observations as follows. If 
we start a trajectory with the same initial conditions in the cases Ah = and Ah small, the latter 
trajectory will shadow (follow approximately) the former for a certain length of time. However, the 
latter will gradually (linearly in time) deviate from the first trajectory due to the (weak) defocusing 
effect of the boundaries, eventually becoming effectively decorrelated. For smaller values of Ah, 
the shadowing will persist for a longer period. 

We refer again to the CTRW model described in Sec. 15.51 For Ah = 0.01, we find that the 
density function Y{t) oscillatory, so that it is difficult to determine its decay rate; instead, 
Fig. l5.27l shows for long times, where 

oo 

:= I Y{t')dt'. (5.23) 

f 

^[This was done in ISL06I . where a simple scaling form was found allowing a data collapse of the data for different 
M.] 
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Figure 5.26: Log-log plot of the diffusion coefficient as a function of A/i for the zigzag model. 

It decays like ~ f"^-^^, so that ~ t^^^^-^^, giving v = 2.23, which according to (15.17b 
gives normal diffusion, as required. However, we may not have managed to attain the asymptotic 
regime in this calculation. For smaller Ah we expect the distribution to have a progressively longer 
tail, so that the average laminar length will be longer, resulting in a larger diffusion coefficient as 
seen in the numerical experiments. 

Again we remark that although the same behaviour is found in the polygonal Lorentz model 
at the level of the statistical properties, it is less clear how to obtain a qualitative understanding in 
that case without identifying the type of laminar behaviour which is presumably responsible. 
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Figure 5.27: Tail of ^{t) on a log-log plot. 



Chapter 6 



Three-dimensional periodic Lorentz gases 



We would like to extend the detailed knowledge available for 2D systems to study more physically 
relevant 3D periodic Lorentz gases. Relatively little work has been done on higher-dimensional 
models, since computation time increases and technical difficulties arise, some of which are dis- 
cussed below. In this chapter we succeed in addressing some of these issues and point out where 
more work is still requirecj^ 

Previous work on 3D periodic Lorentz gases in a physical context includes that of Bouchaud 
& Le Doussal IIBLD85II . who studied the Lyapunov exponent and the decay of velocity auto- 
correlation functions for simple (hyper-)cubic lattices in <i ^ 7 dimensions, and Dettmann et al. 
IIDMR95I . who studied a 3D hexagonal close-packed lattice with an electric field and a thermostat 
(a device for keeping the kinetic energy of a particle constant, despite the energy input due to the 
electric field), showing that the Lyapunov exponents satisfied a conjugate pairing rule. Rigorous 
results, discussed in more detail below, are proved in IIChe94ll and IIGWOOII . 

We first review results which show that higher-dimensional Lorentz gases with finite horizon 
exist. We then construct a particular 3D model with overlapping scatterers which has a finite 
horizon regime, analogously to the model studied in Chap. [3l in this regime we show numerically 
that our model is diffusive. We then consider the effect on the statistical properties of allowing an 
infinite horizon, and find that there are two qualitatively different types of infinite horizon regime. 

6.1. Existence of higher-dimensional Lorentz gases with finite horizon 

Chernov IIChe94ll extended the results of IIBS811 IBSC91II to higher-dimensional {d ^ 3) periodic 

Lorentz gases. He provej^fast (at least stretched exponential) decay of correlations for the billiard 

map, and the central limit theorem and functional central limit theorem for Holder continuous 

observables for the billiard flow, when the scatterers are disjoint and have sufficiently smooth (C^) 

^[Some results from this chapter are included in a submitted manuscript ISanl . which also contains updated refer- 
ences.] 

^Some details of the proofs in IChe94l are corrected in IBCST02IIBCST03I . 
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boundaries, and the finite horizon condition is satisfiecj^. However, no explicit example of a model 
satisfying the assumptions was given in ||Che94i . even in the 3D case. In fact we are not aware of 
any examples of such models in the literature; unfortunately we have also been unable to construct 
an exphcit example. 

6.1.1. Rigorous results from convex geometry 

Several results in convex geometry, which are seemingly unknown in the physics community, 
bear light on the possibility of constructing higher-dimensional periodic Lorentz gases with finite 
horizon. 



In | Hep61[ it was shown that any lattice packing of spheres in 3D has an infinite horizon. 



and in fact a cylindrical hole, the three-dimensional version of the corridors in infinite horizon 
2D periodic Lorentz gases. Such a hole consists, in billiard language, of a collection of parallel 
trajectories (forming a cylinder), none of which ever collides with a scatterer. 

We remark that the term 'lattice' in this result refers to a set of points of the form ^"^j a,e,, 
where a, € Z; thus geometers' lattices are what physicists call Bravais lattices IIAM76II . For 
example, the hexagonal close-packed structure is not a lattice in this sense, although the face- 
centred cubic is. A 'sphere packing' is (roughly) a collection of touching, but non-overlapping, 
spheres. 

The result was later extended to show the existence of cylindrical holes in lattice sphere pack- 
ings of any dimension: see references in MHZOOI . This implies that to obtain a finite horizon it is 
necessary to have more than one sphere per unit cell; this was stated without proof in BCDOOII . 

Recently it was proved in IIHZOOII that this is no longer true if we consider arbitrary periodic 
structures, rather than just lattices, consisting of identical convex bodies. They showed that in any 
dimension n there exist periodic arrays of non-touching spheres with finite horizon (although they 
did not use this terminology), and in fact with spheres replaced by any convex body. The proof is 
constructive, but it is not easy to convert it into the construction of an explicit examplqj, since it 
involves finding sets of minimal cardinality satisfying certain properties. 

This result impUes that the class of models considered in ||Che941l is indeed non-empty. We 
remark that in any such finite horizon model, the minimum number of spheres which can be seen 



Recall from Sec. l3.1l that a periodic Lorentz gas satisfies the finite horizon condition if no trajectory can be extended 
arbitrarily far without hitting a scatterer. 
^M. Henk, private communication. 
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from a given sphere has been shown to be at least 30: for references on this and related results see 
the review llZonOll Sec. 5-6]. Hence any such structure must be quite comphcated. 

6.1.2. Attempts to construct a finite-horizon model 

A common construction in solid state physics is to build lattices as stacks of layers, with each 
layer being a 2D triangular lattice of spheres (corresponding to the standard triangular 2D periodic 
Lorentz gas). Consecutive layers are placed in one of three possible positions, labelled A, B and C, 
relative to the previous two layers IIAM76II : a hexagonal close-packed structure coiTcsponds to the 
choice . . .ABABAB . . ., and the face-centred cubic to . . . ABC ABC . . .; neither of these (somewhat 
surprisingly) has a finite horizon. (The face-centred cubic structure is a Bravais lattice, and so is 
covered by the theorem cited above, but the hexagonal close-packed structure is not.) We could 
ask if it is possible to produce a finite horizon by some other periodic ordering of As, Bs and Cs. 
However, in fact placing adjacent layers on top of each other creates cylindrical holes between 
the layers, so that there is always an infinite horizon. If we also require disjoint (non-touching) 
scatterers, then we would need to separate the scatterers slightly from their touching positions, 
thereby introducing more holes. 

As an example. Fig. 16.11 shows the case of the face-centred cubic lattice: Fig. 16.1 f a) shows 
the structure as built up from stacks of hexagonal close-packed layers, whilst Fig. I6.1f b) shows 
a cylindrical hole through the structure, as required by Heppes' theorem quoted above. In fact 
this hole lies between two adjacent hexagonal layers, so that any structure built up of such layers 
would contain a similar cylindrical hole. 

Another possibihty for constructing explicit 3D models with finite horizon would be to make 
periodic layers with a lower density of spheres, but still with finite horizon within the layers. 
If the density were low enough, the next layer could be pushed down far enough to block the 
corridors between the layers and possibly obtain a finite horizon. However, we have not managed 
to construct a model in this way. We could also try constructing a model with, for example, 
ellipsoids instead of spheres, but again we have had no success. 

6.2. Construction of 3D periodic Lorentz gas with overlapping scatterers 

Here we present a 3D periodic Lorentz gas model analogous to the 2D modejf] of Chap. [3l which 
has a finite horizon in certain regions of parameter space. To accomplish this, however, we allow 
^In fact, we were led to the 2D model as a cross-section of the 3D model. 
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(a) (b) 

Figure 6.1: (a) Face-centred cubic lattice as constructed from hexagonal close-packed layers of touching 
(but non-overlapping) spheres, (b) Two small cylindrical holes (top left) between two layers. 



the spheres to overlap, which appears to be the only way of achieving a finite horizon with a 
simple model, as discussed above. 

Allowing the spheres to overlap is at first sight non-physical; however, suppose that the moving 
particles are spheres of non-zero radius s. In this case, the dynamics is the same as that of point 
particles moving through the same lattice, but with the radius of the scatterers increased by s; this 
radius could represent the effective radius of a particle interacting with the crystal lattice via a 
more realistic potential. This is the same construction originally used by Sinai IISin70ll to reduce 
two discs moving on a torus to a periodic Lorentz gas; see also the discussion in Chap.[T] We will 
find conditions under which the overlapping model can be regarded as a physical model in this 
way. 

Note that allowing overlapping scatterers violates one of the key requirements in Chernov's 
proof that the multi-dimensional Lorentz gas is diffusive IIChe94ll . Nonetheless, we expect phys- 
ically that the system still possesses strong ergodic and statistical properties, and in particular is 
still diffusive. This seems to be corroborated by numerical experiments, as discussed below. 



6.2. Construction of 3D periodic Lorentz gas with overlapping scatterers 



139 



6.2.1. Construction of finite iiorizon model by blocking corridors 

We consider a cubic unit cell of a 3D simple cubic lattice of spheres of radius a. If the spheres do 
not touch, then there are lots of holes in the structure. In order to construct a finite horizon model, 
we need to block all of these holes. We first allow the a spheres to overlap. The configuration in 
a plane through the centres of 4 neighbouring a-spheres then looks like Fig. |6.2(a)[ There is no 
longer an infinite horizon within this plane, and we have one corridor perpendicular to the plane. 

We now add a new scatterer, with (different) radius b, at the centre of the a-unit cell. If we 
make the Z^-sphere sufficiently large that its projection onto one of the faces of the unit cell blocks 
the hole in that face (Fig. |6.2(a)l ), then we have blocked all corridors perpendicular to the faces. 

6.2.2. Phase diagram 

The construction of the phase diagram follows that of the 2D case: again we look for lines in pa- 
rameter space across which the geometry of the unit cell undergoes qualitative changes. However, 
the additional freedom in 3 dimensions makes the calculation more difficult. 

Overlap of a-spheres The overlaps of the a-spheres change as follows when we vary a. 

• The a-spheres overlap if and only if a ^ ^ ; this is a necessary condition for a finite horizon 
in our class of models. 

• The space between the a-discs on a unit face closes when two discs at opposite ends of 
a diagonal of the square touch, which occurs when = for larger values of a it is 
impossible for a particle to move between different unit cells, so the trajectory is locahsed 
and the diffusion coefficient vanishes. 

• The a-spheres cover the entire unit cell when two spheres at diametrically opposite vertices 
of the unit cube touch, which occurs when 2a = \/3, i.e. at a = 

Blocking vertical trajectories Let bmin be the minimal value of the radius b of the additional 
scatterer required to block the vertical corridor (i.e. such that its vertical projection covers the space 
on a face between the a-sphere overlaps, as in Fig. |6.2(a)| ), and let d be the width of the overlap 
of two a-discs on a face; this is also the radius of the disc of intersection of two neighbouring 
a-spheres with the mid-plane between them. The geometry is shown in figure [6^(a)| 
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(a) 

Figure 6.2: (a) Geometry of overlapping discs on one (and hence each) face of the cubic unit cell. The 
dashed circle is the projection of the ^-disc at the centre of the cell, (b) Geometry in the mid-plane of a unit 
cell. 

I and a necessary condition for 

(6.1) 

Blocking diagonal trajectories We now must now block any diagonal corridors (diagonal rel- 
ative to the lattice viewed as cubic cells). By symmetry, it suffices to block trajectories in the 
mid-plane of the unit cell: there is 'most space' in this plane. (Compare this to blocking diagonal 
angles at 45° in 2D.) 

The mid-plane is shown in Fig. |6.2(b)[ its geometry is the same as that of the 2D model 
considered in Chap. [3j after replacing aby d. If b ^ ^min> then horizontal and vertical trajectories 
within this plane are already blocked, by definition of ^min- Diagonal trajectories at 45° will be 
blocked if either d ^ which reduces to a ^ 0.612, or if ^ 

If neither of these conditions holds, then, by continuity, there is an interval of non-zero size 
around the height of the midplane such that planes at heights lying in that interval also have an 
infinite horizon along trajectories at 45°, so that there is a cylindrical hole of non-zero volume 
lying along this diagonal, and hence the structure has an infinite horizon. 



We have Z^min + d = ^ and + (^)^ = a^, so that d = \J ■ 
finite horizon is 



6.2. Construction of 3D periodic Lorentz gas with overlapping scatterers 



141 



Conditions for localised trajectories If the a or b spheres are too large then they will overlap 
to such an extent that trajectories will be localised, as in the S and D regimes of the 2D model 
(Sec-ETl- 

Consider again the midplane shown in Fig. |6.2(b)[ Suppose that the b- and (i-discs do not 
touch. Then any point outside the discs in the midplane is accessible from any other. In fact, since 
there are identical planes perpendicular to each of the coordinate directions, there is a connected 
network of empty space lying around the grid of Unes parallel to the coordinate directions which 
join the centres of the ^-spheres. Hence for b + d < we do not have locahsed trajectories. We 
now show that in fact this is a necessary condition: all trajectories are localised if b ^ -^ — d = 

Consider cutting the unit cell by a horizontal plane which begins at the mid-plane and moves 
upwards. Let h be the height of the plane above the mid-plane, and let the radius at height h of the 
a- and Z7-cross-sections be a{h) and b{h), respectively. Then a(0) = d, b{0) = b, b{h) = ^Jb^ — h^ 
and a{h) = \j — {\ — ^Y- 

Trajectories will be localised if and only if there is a plane at some height h in which there is 
no available space, i.e. for which the aih)- and Z7(/i)-discs fill this plane, for then there are barriers 
in each direction blocking escape, whilst if there is space in each plane h then in fact this space 
must be connected all the way from /j = to /j = From Sec. 13.1.41 we know that this blocking 
occurs exactly when 

b[h) ^ \ - ^Ja{hY-l (6.2) 

since the geometry in this plane is the same as the 2D geometry with a and b replaced by a{h) 
and b{h), respectively. Considering when there exists an h for which b{h) = ^ — \J a{h)^ — ^ is 
equivalent to the following slightly different approach. 

Taking coordinates with the origin at the centre of the unit cell (and hence at the centre of the 
Z?-sphere), if the intersection point of the overlap of neighbouring a-spheres with the ^-sphere is at 
{x,y,z) then we have ;ic = by symmetry and thus 

+ = b^, and (6.3) 
(0-i)2 + (3;_i)2 + (,_i)2 = ^2^ (6.4) 

since the intersection point lies on each sphere. 
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Expanding (16.4b and using (16.31) gives y + z = cc '■= b — a + ^, and substituting z = cc—y into 
631) then gives 

y = Ua± - a^] , (6.5) 



2 

so that an intersection exists if and only if 2b^ — ^ 0, which is equivalent to b ^ since 
b>0. 

We have t? < so that < |, and hence a = b^ — {a^ — > 0. Thus the condition for 
existence of the intersection is b\/2 ^ b^ — + ^, i.e. b^ — bV2 + (| — a^) ^0. Equality occurs 
when b = ^± ^Jd^ — \, so that there is an intersection if and only if 



The leftmost term is equal to — d; the rightmost term is greater than so that it is never 
attained for non-localised trajectories for which we must have a < Hence the condition for 
intersection is that the b- and J-discs touch in the mid-plane, as claimed. 

Condition for all space to be covered To establish when the Z7-sphere covers all the space left 
by the a-spheres, consider first the case when a < so that the a-spheres leave holes on and 
near the unit cell faces. These holes will be covered by the Z7-sphere once it is big enough that its 
intersection with the unit cell face covers the space on that face, i.e. if 

b{h)^b,^^ for/j = i (6.7) 



which reduces to 



Z,2-(l)2^ 1 /«2_1 (6.8) 



or finally 



(6.9) 



When a > we need b to be large enough for the Z^-sphere to cover the 3-hole inside the 
overlap of the a-spheres. For this we need that the Z^-disc covers the 2-hole between the J-discs in 
the mid-plane. The geometry is exactly as in the 2D case, with a replaced by d, so we can quote 
the result: 



^i- V^2_l = 1_ /^2_l. (6.10) 
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Overlapping conditions Note that ifb> \/3/2 — a then the Z^-spheres meet the a-spheres, whilst 
if b > 1/2 then neighbouring Z?-spheres meet each other. 

Phase diagram The parameter space of the 3D model is shown in Fig. 16.31 for b < a. The 
regimes are named as far as possible to agree with those of the 2D model in Sec. 13.1.41 if the name 
is the same as one for the 2D model, then the geometrical features and statistical properties are 
similar. The diagonal from top left to bottom right is the boundary line separating regimes where 
the Z7-spheres do not touch the a-spheres (below) from those where the b- and a-spheres overlap 
(above). This cuts several of the other regimes. 

Several regimes do not occur in the 2D model. N has localised motion in cells centred on the 
Zj-spheres, but with the Z?-spheres overlapping the a-spheres. In IH4 the a-spheres meet each other 
and also the Z7-spheres but still with infinite horizon. In IHO nothing is touching and we have a 
strongly infinite horizon (see also Sec. 16.51) . 

Note that the finite horizon (FH) regime consists of two disjoint pieces, corresponding to the 
two different ways of preventing diagonal infinite horizon trajectories detailed above. In fact, the 
top of FH2, above the line b = could be regarded as a third finite horizon regime. Here, the 
Z^-spheres meet each other (as well as the a-spheres), so that there is a single connected scatterer 
network with an interconnected hole inside it. 

Figure |63] shows the model in the finite horizon regime with the a- and Z7-spheres not overlap- 
ping. We can check by rotating the model that it does indeed have finite horizon, since no holes 
can be seen in the complete unit cell. 

6.2.3. Physical realisation with a moving particle of non-zero size 

The overlapping system with radii a and b studied above can be regarded as physical if it is a 
reduction of a system with a moving sphere of radius s travelling within a lattice of spheres with 
radii a' and b' with no overlaps. Regarding {a,b) and {a',b') as position vectors in the phase 
diagram determined above, the reduction condition is 

{a,b') + {s,s) = {a,b), (6.11) 

which corresponds geometrically to translating by the vector {s,s) from the initial point on the 
phase diagram. The system will be diffusive if this lands in the FH regime. 
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Figure 6.3: Parameter space of the 3D model showing the different regimes described in the text. 

For the initial point to correspond to a system where no spheres overlap, we must start in the 
region IHO. In fact, any point in IHO corresponds to some point in FH via a translation by some 
vector {s,s), except for a line at 45° passing through the point where the FHl and FH2 regions 
touch. 

6.3. Normal diffusion in the finite horizon regime 

We now use techniques developed in previous chapters to give strong evidence that we have nor- 
mal diffusion in the finite horizon regime. We emphasise that the arguments of ||Che94|| require 
disjoint scatterers, and so do not immediately apply to our model; nonetheless we expect that those 
methods should be able to be refined to deal with our case and prove normal diffusion, even in the 
strongest sense that the functional central limit theorem is satisfied. 

6.3.1. Decay of velocity autocorrelations 

From the discussion of the Green-Kubo relation in continuous time of Sec. 12.3.81 a necessary 
condition for the existence of the diffusion coefficient (which is the weakest type of normal diffu- 
sion) is the integrability of the velocity autocorrelation function C{t). Since this function is highly 



6.3. Normal diffusion in the finite horizon regime 



145 




Figure 6.4: The 3D Lorentz gas with finite horizon in a regime (r = 1.56, b = 0.31) where the b- and a- 
spheres do not touch, as can be seen in (a); a different view of the structure is shown in (b), which confirms 
that there is a finite horizon, since no holes can be seen in the complete unit cell. 
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oscillatory and rapidly becomes small, it is difficult to study. It was thus suggested in IILM93II to 
look instead at the integrated velocity autocorrelation function R{t), given by 





The integration in the definition of R{t) acts as a smoothing operation which deals efficiently with 
the highly oscillatory nature of C{t). Note that in this definition we do not need to perform a direct 
numerical integration of the function C{t); such an approach was taken, for example, in [|MM97II . 

If C{t) is integrable, then R{t) dD as ? ^ oo, where d is the spatial dimension and D is 
the diffusion coefficient. We find numerically that C{t) decays so fast that we cannot extract any 
information about it, and we do not plot it. Similarly R{t) reaches very rapidly a limiting value, 
which it then oscillates around in a random fashion. This provides evidence that C{t) is integrable 
and hence that D exists, giving normal diffusion. 

6.3.2. Growth of moments 

Figure [63] shows the growth rates 7^ of the moments (|r|^) as a function of q. As in the 2D 
finite horizon case (not shown), they satisfy 7^ = qjl, so that we have a strong form of normal 
diffusion; we conjecture that the functional central limit theorem is satisfied for our model in the 
finite horizon regime, as was proved in ||Che94|| in the case of disjoint scatterers. 

6.3.3. Shape of distributions and central limit theorem 

Knowing that the mean squared displacement grows asymptotically linearly in the finite horizon 
regime, we can apply the techniques of Chap.|4]to investigate the shape of densities and the central 
limit theorem. 

Restricting attention to ID projections of the position distribution, we need to calculate the 
analogue for our model in the finite horizon regime of the available height function h(x) used 
in Chap. IH This analogue is the available area, which we denote by A(x), in a cross-section 
perpendicular to the x-axis at distance x from the centre of the Z7-sphere; we refer to this cross- 
section plane as the x-plane. We denote the radii of the cross-sections of the a- and Z^-spheres in 
that plane by a(x) := ^ — {x—^Y and b{x) := Vb^ —x^, respectively, where we again use the 
convention that ^/a = when a < 0. 




(6.12) 
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Figure 6.5: Growth rate 7^ of the moments {\r\'') for particular parameters in the finite horizon regime 
ir = \ .6\ b — 0.3), as a function of q. The growth rates obey well the relation Jq = q/2, ?& for the 2D 
Lorentz gas with finite horizon. 



Figure |6.2(a)| can be thought of as depicting a generic cross-section, after possibly adding a 
small Z?-cross-section in the centre. Denote by Aoveriap (■'c) the area of each overlap of a(x)-discs in 
the jc-plane. Then 

A{x) = 1^ [1 - T^H^f+K^f) +4AoverlapW] , (6.13) 

where if a(x) > 1/2 then 

AoveriapW= j 2^ a{xY - dy = 2[^a{xf - F {{-x)] , (6.14) 

y=l/2 

whilst there is no overlap ifa{x) ^ 1/2. Here 

F{y;x) -yJaixJ^-y^ + ^aixftm-' ( —=!==] (6.15) 

Vv«w -r/ 

is the anti-derivative of \J a{x)^ —y'^ (with respect to y). Further, |Q| is the volume of the available 
space in a unit cell given by 



el = l-|7r(a^ + Z73)+3Voverlap, 



(6.16) 
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Figure 6.6: ID position density ft{x) at time f = 50 for the 3D Lorentz gas in the finite horizon regime 
at r = 1.61, = 0.3 (for which the a- and ^j-spheres do not overlap). The inset shows the available space 
function A (jc), and the main figure also shows the demodulated density Pt{x), compared to a Gaussian with 
variance 2Dt, where D ~ 0.081 for this geometry. 



where 



Voveriap:=2 J K{a^ - j^) Ay = Ik [la^ - kci^ + M 



(6.17) 



1/2 



is the volume of the intersection of two neighbouring spheres. (Similar overlap calculations were 
used in MSanOOII in a related context.) In the above calculations we have for simplicity restricted 
the calculation to the case where the b- and a-spheres do not overlap. 



The inset of Fig. [63] shows the available space function A{x) as calculated above, whilst the 
main part of that figure shows the ID position density ft{x) and the demodulated version f\t{x) := 
ft{x)/A{x), for certain geometrical parameters in the finite horizon regime which we believe to be 
representative. Again we see that the demodulated density is very close to Gaussian, and we have 
confirmed numerically that the central limit theorem holds with rate of convergence close 
to the optimal t^^l'^ as discussed in Sec. 14.61 see Fig. 16.71 



6.4. Statistical behaviour in the infinite horizon regime 



149 




^'1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 



Figure 6.7: Convergence to limiting normal distribution with variance 2D for r = 1.61, b = 0.3. The rate 
for long times is f^'*'^^. 



6.3.4. Geometry dependence of diffusion coefficients 

Figure [Q] shows the geometry dependence of the diffusion coefficient over the two finite horizon 
regimes. In FHl the behaviour resembles that in the 2D case, although the angle at which the 
curves approach on the right hand side is different; in particular, here too there is a qualitative 
change in behaviour with a non-trivial maximum which disappears for larger values of r. 



6.4. Statistical behaviour in the infinite horizon regime 
6.4.1. Shape of 2D distributions 

How can normal diffusion fail to hold in dispersing billiards with infinite horizon? Figure 
shows scatterplots representing the position and displacement distributions for representative pa- 
rameters in the infinite horizon regime of a 2D square Lorentz gas. These distributions have a 
distinctive shape caused by the possibility of particles travelling arbitrarily far without ever hitting 
a scatterer. Bleher ||Ble92|| showed, assuming some natural conjectureCl that the mean squared 



^[A complete, rigorous proof has now been given ISV07I .1 




Figure 6.8: Diffusion coefficients for parameter values in the finite horizon regimes: (a) in FHl; (b) in 
FH2. In (a), r increases from r = 1 .425 to r = 1 .625 from bottom to top, in steps of 0.025. In (b), r increases 
from r = 1 .65 to r = 1 .975 from bottom to top, in steps of 0.025. 
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displacement in this case grows like (Ax^)f ~ tlogt as f ^ oo. He also showed that 

Xf-xo V 



z, (6.18) 



where z is a normal random variable; this is a central limit theorem type result with a different 
normahsation constant. The faster growth rate of the variance corresponds to the tails of the 
distribution visible in Fig. 16.91 

6.4.2. Discrete-time dynamics 

Following ||Ble92|| (who discussed the 2D case), we write 

«— 1 n— 1 

X„-Xo = 52(Xy+l -X;) = Xl^'i' 

;=o ;=0 

where x„ := x(f„) is the position at the nth collision (occurring at time tn) and := xy+i — xj is 
the free flight vector between the jth and [j + l)th collisions. 

In the finite horizon regime the free path length is bounded above, so that is a piecewise 
Holder continuous function defined on the collision phase space. The central Umit theorem of 
IIBS811 IBSC91II can then be applied to show that (x,, — xo)/a/« converges in distribution to a 
normal distribution. 

In the infinite horizon case, however, rj is no longer bounded and so cannot be piecewise 
Holder continuous, since a continuous function on the compact phase space M must be bounded. 
The central limit theorem thus does not apply. Defining the finite-time diffusion coefficient D„ by 

D„:=-^(|x„-xo|^)v, (6.20) 

the limiting diffusion coefficient for the discrete time dynamics takes the form ||Ble92|| 

1 1 °° 

D^:= limD„ = -(|ro|') + -£(roTy), (6.21) 

if it exists. 

In ||Ble921 it was conjectured that (at least in 2D) (ro • ry) decays fast for large j, so that 
the infinite sum in (16.211 ) exists. However, it was proved there that (|ro|^)v = ('7^)v, the second 




Figure 6.9: (a) 2D position distiibution; (b) 2D displacement distribution, r = 2.5; b = 0.0; / = 50; 
N = 5 X 10^ initial conditions. 
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moment of the free path length for the billiard map, diverges (logarithmically), so that the first term 
in (16.211 ) is infinite and hence so are the discrete-time and continuous-time diffusion coefficients. 

6.4.3. Obstructions to normal diffusion 

Generalising the above to higher dimensions, there are two possible obstructions to normal diffu- 
sion outside the FH regime: 

(i) the continuous-time VACF C{t) decays slowly as a function of time t, so that the Green- 
Kubo integral, giving D as the infinite-time integral of this VACF, diverges; and 

(ii) the second moment of the (discrete-time) free path length diverges. 

If neither of these obstructions is present, then we expect to have normal diffusion: this is certainly 
the case in the finite horizon regime treated above. We now discuss the situation in the infinite 
horizon regime, i.e. when there exist trajectories which never collide with the scatterers. 

6.4.4. Review of two-dimensional case 

Decay of velocity autocorrelations The following type of argument seems to have first been 
published in IIFM84I . Consider a 2D lattice with infinite horizon (for example a square lattice of 
scatterers of radius a). Due to the infinite horizon, there exist corridors in the structure, i.e. empty 
regions of constant, non-zero width, extending infinitely far in opposite directions. Bleher ||Ble92|| 
gives a detailed description of these corridors. 

Consider one particular corridor, of width d > 0. Take coordinates such that the y-axis is 
parallel to the centre-line of the corridor and the x-axis is perpendicular to it, with the origin at an 
arbitrary point on one edge of the corridor. Let Q be the angle from the y-axis. 

We consider initial conditions (xo,3'o) in whose trajectories have no collisions within (con- 
tinuous) time t, and in fact remain inside the corridor at least until time t; thus these trajectories 
are straight during this period, so that Vj = Vq. A trajectory emanating from (xcjo) in the direction 
Q from the j-axis will not leave the corridor within time t provided the line segment of length t 
at angle Q ends within the corridor, i.e. provided 0_ ^ ^ 0+, considering only one direction of 
flight down the corridor initially. The geometry gives 



sin0 



d — x 



sm d = — 



(6.22) 



'+ = 



t 



t ' 
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The Green-Kubo formula for the diffusion coefficient means that we are interested in the rate 
of decay of the velocity autocorrelation function (VACF) ast ^o°. For large t, 1/f is small, hence 
sin 0-1- is small, so that sm6± 6±. 

For large enough t, denoting the average of the VACF solely over those long trajectories which 
do not escape after time t by (vq • Vf)iong, we thus have 



Here, K and K' are constants which are related to the area of intersection of the corridor with V. 

It is then argued in IIFM84I that if the rest of the trajectories are more efficiently mixed, as 
we expect they are due to collisions with the scatterers, then the VACF averaged over those other 
trajectories will decay faster. Hence the rate of decay averaged over the whole of V will be 
dominated by the slow l/t decay of the long trajectories considered above. In this case, the 
integral of the VACF will diverge, so that the diffusion coefficient will not exist (or is infinite), 
implying that the system is super-diffusive. 

Unfortunately the above argument is not rigorous, and indeed currently there are no rigorous 



results on decay of correlations for the billiard flo\M3, even in the finite-horizon case HCYOOI . 
although the recent results of IISV03II prove that £) = oo, so that they show indirectly that C{t) must 
be non-integrable. 

In IIGG94II the exact result of Bleher for the rate of increase of the mean squared displacement 
was compared to numerical results. They pointed out, and we have confirmed, that the coefficient 
of the flog? term is very small, so that it is difficult to observe numerically. Nonetheless it is 
sometimes possible to observe this effect by plotting {Ax^)t/t against logf, as first suggested in 
the context of ID maps in IIGT84II . As argued in ||GG94i if {Ax^)t ~ flogf then the lower order 
terms presumably take the form 




d 



{d-x)/t 



(6.23) 



x=0 



e=-x/t 




(Ax^), r^A + Bt + Ctlogt. 



(6.24) 



It then follows that 



Ar^ +B + Clogt. 



(6.25) 



[Such results have now been obtained in the finite-horizon case 
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Introducting the variable z '■= logf, we then have 

^-Jl^Ac-''+B + Cz, (6.26) 

so that for large enough values of z, i.e. for large enough values of t, we should find asymptotically 
linear growth of {Ax^)t/t. 



Growth of higher-order moments In IIAHQ03L a similar line of argument was used to give 
a lower bound on the rate of increase of higher order moments, as follows. A proportion C/t of 
trajectories do not collide in time t (as described above), so that for those trajectories we have the 
lower bound 

{\rniono^jivtr = Ct''-\ (6.27) 
where v is the speed, so that 7^ ^ ^ — 1. 

They show that 7^ is a convex function of q. We also know that 7o = and 72 = 1 lie on the 
curve, that ^ q/2, assuming that the process acts at least as fast as a random walk, and that 
Yq ^ q- These together lead to the conclusion that 



q/2, when q <2 

(6.28) 

q—l, when ^ > 2. 



Tail of free path distribution As discussed in Sec. 16.41 the existence of Ey [x^] is crucial for 
the possibility of having normal diffusion. We introduce the following functions: 

»p(r) := Pv (t > r) := V {{x G M: r{x) > T}) ; (6.29) 
»F(r) := (t >T):=n{{xeM: z{x) > T}) . (6.30) 

(Recall that the bilhard map T preserves the measure v on the phase space M, and the billiard 
flow preserves the measure jJ. on the phase space Ai.) These functions describe the tail of the 
distribution of the free path length considered in discrete time and continuous time, respectively. 
Note that in the discrete time case |r| = t, where r is the free flight vector. 
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In 2D, it was proved in ||Ble92|| that 

»F(r)~r"2 asT^oo. (6.31) 

By definition of expectation, 

r 

Ev [t^] = - lim [ d»F(r), (6.32) 



where the integral is a Lebesgue-Stieltjes integral. Integration by parts is valid for such integrals, 
so that 

- jT^d'¥{T)= j2T'¥{T)dT + 0{\)r^logT', (6.33) 

7b To 

and hence the expectation diverges logarithmically. 

The same result can be obtained using the following relation proved in IIDDG97II (see also 

CvEv[/(t)]=c^E^[/(t)], (6.34) 

for any function / : R+ —>■ M such that /(O) = 0; here /' denotes the derivative of /. For 
f{z) = z we recover the known expression (see IIChe97ll ) for the mean free path in terms of Cy and 
discussed in Chap.[3l while for f{z) = we obtain 

[t] = ^Ey [z^] . (6.35) 

Here we extend the definition of t to the whole of Ai as 

t(x,v) :=inf{f gM+: x + fveM}, (6.36) 

the minimum time needed to collide with a scatterer. 



Thus the second moment of T with respect to v exists if and only if the mean of T with respect 
to jj. exists. We now use an argument similar to those above involving the proportion of trajectories 
which have not collided within time t: those trajectories have (continuous-time) free paths which 
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are at least t, so that integrating over initial conditions in one unit cell we have 

'i'{T)^CT \ (6.37) 
The mean [t] thus diverges logarithmically, and hence so does Ey [t^] . 

Summary The above arguments all involve calculating the proportion of trajectories starting 
from a single unit cell which undergo no collisions up to time t. For large t these trajectories are 
ones which remain inside the corridor and lie on a circle of radius t. 

We obtain rigorous lower bounds for the higher-order moments (|r(f)|^) as f ^ oo and for the 
size of the tail *F(r) of the free path distribution function as T ^ oo, whilst we obtain only a 
heuristic estimate of C{t) as f ^ oo. In the next sections we shall apply these methods to our 3D 
Lorentz gas. 

6.5. Simple cubic lattices 
6.5.1. Free path distribution 

We begin by discussing simple cubic lattices with r > 2 and b = 0. In IIChe94ll it was stated 
that point (ii) of the previous section, i.e. divergence of Ey [t^] holds in higher dimensions, and 
that the proof is "easy to verify", but it was not stated which configurations this applies to. For 
simple cubic lattices a proof was given in HGWOOi In the next section we show that in certain 
other geometries the result is in fact no longer true, so that the existence or otherwise of normal 
diffusion reduces in that case to determining if (i) holds. 

The method of BGWOOI to prove this result is to generalise the method discussed in the previous 
section, by using the existence of free planes (called 'sandwich layers' in liGWOOl ) in the simple 
cubic lattice. These are infinite planes which do not intersect any scatterer, and are the 3D analogue 
of corridors related to infinite horizon trajectories (which we can think of as 'free lines') in 2E(^. 

Analytical argument Analogously to the calculations in the previous section, we need to find 

the proportion of trajectories remaining inside a sandwich layer and lying on a sphere of radius 

t. The intersection of these two objects is approximately a circle with non-zero width, and the 

proportion is then the ratio of the area of this thickened circle, which is IntAx, divided by the area 

^[Figures showing the different types of holes in 3D, as well as an analytical argument for systems of any dimension, 
are given in a submitted manuscript ISanM 
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Figure 6.10: Decay of tail of the free path distribution function for the billiard map for simple cubic lattices 
(b — 0) with r = 3.0 and r = 2.3, compared to the predicted decay rate of t^^. 



4nt^ of the sphere of radius t. The ratio is hence C/t. A more careful argument, taking account of 
the intersections of such circles corresponding to different sandwich layers, shows that the same 
result holds BGWOOII . We remark that this result was also stated in IIFM84II . 

The argument of the previous section now gives a lower bound on the size of the tail *P(r) of 
the free path distribution for the billiard flow of so that again the expectation [t] diverges. 

Numerical calculation Figure l6?T0l shows the tail *P(r) of the free path distribution for the 
billiard map, for two simple cubic lattices. This was obtained by recording only values of the free 
path exceeding a certain lower threshold, since otherwise the predominance of small values hides 
the information about the tail. A straight line corresponding to a decay rate of T^^ is also shown. 
Although the data appear to decay slightly faster than this, the decay rate is decreasing for large 
T; we believe that the asymptotic behaviour does obey the analytical prediction. 

6.5.2. Higher-order moments 

The growth rate Yq of the higher-order moments is shown in Fig. |6.11[ They are in good agreement 
with the analytical prediction, which in this case is the same as in the 2D case reviewed in the 
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previous section: 



q/2, when q <2 
q — when q>2. 



(6.38) 



6.5.3. Decay of velocity autocorrelations 

Figure I6TT2] shows the velocity autocorrelation function (VACF) C{t) as a function of time t. A fit 
to the less-noisy central part of the graph (also shown) gives a power law decay with an exponent 
which is close to 1. In Fig. l6.13l we also show R{t), the integrated velocity autocorrelation function. 
If C{t) ~ as f ^ oo, then we expect that 

R{t)r^logt, (6.39) 

so that Fig. 16. 13l plots R{t) as a function of logf. The hnear growth for long time provides con- 
firmation of the non-integrability of C. We note that C{t) appears to converge to monotonically 
from above. 
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Figure 6.13: Integral R{t) of the velocity autocorrelation function as a function of log? for the simple cubic 
lattice with r = 3.0. 
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6.6. Lattices with cylindrical holes 
6.6.1. Tail of free path distribution 

Analytical argument In the previous section we showed that the existence of free planes in the 
structure impUes the divergence of Ey [t^] , and hence the non-existence of the diffusion coef- 
ficient. However, for our model it is possible to block all free planes, either with the a-spheres 
overlapping, or for large enough values of b with the a-spheres non-touching, leaving only holes 
of cylindrical type. The previous argument proving the divergence of Ey [t^] now no longer holds, 
and in fact this quantity is finite, as follows. 

Again we calculate the proportion of non-colliding trajectories of length t, now which remain 
inside a cylindrical hole of radius r. The set of these allowed directions along the cylinder is a 
circle of area (approximately) nr^, independent of t, whilst the set of all possible directions is the 
surface of a sphere of radius t, with total area Anf'. (Actually there are two such circles in opposite 
directions along the cylinder.) The proportion of such trajectories is thus C/f^, compared to C/f 
in two dimensions and in the case of simple cubic lattices. 

This gives a decay rate of *F(r) ~ T^^ for the distribution function of the free path length for 
the billiard flow, so that E^ [t] exists and hence by (16.351 ) we also have Ey [t^] < °°. Thus point (ii) 
of Sec. l6.4l is no longer an obstruction to the possibility of normal diffusion. This does not appear 
to have been observed previously, although we later discovered that it was stated in ||BGW98t (see 
also IGWOOl ) that there is a maximum decay rate of ^{T) of 1 /T^^^ in d space dimensions. No 
details of the derivation were given, although presumably the argument was based on a similar 
idea. We however give explicit examples of models where this optimum value is attained. 

Numerical results 

[The numerical results in this section are incorrect. This is due to the use of an incorrect method 
for generating the initial velocity distribution. The initial velocities should be generated uniformly 
on a unit sphere (since the speed is fixed to be 1). The incorrect method used for the numerical 
results depicted here was to distribute uniformly the spherical angles d and ^, and then assign 
the components of the velocity based on these angles. This method, however, does not produce 
a uniform distribution on the sphere, but rather produces a noticeable concentration of directions 
close to the poles of the sphere. These concentrations align with the holes in the structure, thus 
skewing the results sufficiently to make it impossible to see the expected effect. 
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Figure 6.14: Tail of free path distributions for the billiard map, for two geometries in which the only holes 
are cylindrical. For r — 1.75, the scatterers overlap, but for r = 2.08 and b — 0.625 they are all disjoint. 
There is good agreement with the analytical prediction of the decay rate T^^. 



Correct numerical results are given in a submitted manuscript MSani There, the initial velocities 
are correctly distributed uniformly on the sphere, by choosing the velocity vector v uniformly in 
[—1, 1]^, rejecting those v with |v| > 1, and then normalising v to unit length. This generates points 
V uniformly on the unit sphere. 

We find that diffusion is indeed normal, as the heuristic argument predicts, if the only holes are 
cylindrical. ] 

Figure \6AA\ shows the tail of the free path distribution for the billiard map for two geometries 
for which the only holes are cylindrical. The tails decay as T^^, so that Ey [t^] < oo, in agreement 
with the analytical calculation. We remark that the step-like character visible in Fig. 16.141 corre- 
sponds to the fact that the mass of the distribution is localised near points of a periodic lattice, as 
remarked for the 2D case in ||Ble92i We do not expect to see this effect in the continuous-time 
free path distribution, and indeed it is absent in plots of that distribution given in liGWOOl . 
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6.6.2. Higher-order moments 



The generalisation to 3D of the argument of IIAHQ031 discussed above gives 



(6.40) 



This does not give enough information to determine uniquely the shape of the curve as it did in 
two dimensions, although it does show that it must grow like q for large q. Numerically we find 
that the data, shown in Fig. 16.151 fit well the function 



q/2, when ^ < 3 
q — Ji/l, wheng>3. 



(6.41) 



The occurrence of the number 3 here is perhaps related to the fact that we are now in <i = 3 space 
dimensions. 
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6.6.3. Velocity autocorrelation function 

The decay of the velocity autocorrelation function at long times has a lower bound of C'?^^, by 
the 3D version of the Friedman-Martin argument discussed above. It is thus now possible that 
the VACF C{t) is integrable and thus that the diffusion coefficient could exist, despite the infinite 
horizon. 

Figures [6. 1 6146. 1 8 1 show the growth of the integrated velocity autocorrelation function R{t) := 
fQC{s)dLS for several sets of geometrical parameters. As argued in the previous section, if C{t) ~ 
t^^ as f ^ oo then we expect R{t) ^ log?. The data in Fig. 16. 16l provide strong evidence that for 
r = 1 .7 and = 0.0, we do have such a f ^ ^ decay, and hence that the diffusion coefficient does not 
exist in this case. 

Figure lOT] shows R{t) for r = 1.8 and b = 0.3. This differs from the previous plot in the 
presence of a large central scatterer which severely restricts the size of the cylindrical holes. The 
figure is not conclusive, but it is possible to believe that R{t) converges as f ^ oo, implying that 

D < oo. 

The most interesting case is that for r = 2.07, b = 0.6, shown in Fig. 16.181 since here the 
scatterers are disjoint. The data here seem to be consistent with R{t) ~ \ogt, and hence C{t) ~ t^^, 
although perhaps the growth levels off for larger values of t. We have also looked directly at C{t), 
but again this converges rapidly to zero. There is some evidence that it oscillates around zero, 
which gives more chance of having convergence. 

These remarks suggest that there may be a possibility of having normal diffusion, in the sense 
of asymptotically linear growth of the mean squared displacement, even in the presence of an 
infinite horizon, provided that there are only cylindrical holes; and this may even be possible with 
disjoint scatterers. 

However, our numerical evidence also points to the fact that many parameters for which there 
are 'large' cylindrical holes in fact have C{t) ~ resulting in super-diffusion. (We remark 
that extending our numerical results to longer times would be difficult, since the calculations 
for Fig. 16. 181 required a total of approximately 6 weeks of CPU time on modern workstations, 
distributed over 16 processors.) This slow decay is not predicted by the heuristic arguments given 
above, and it is important to determine its origin. One possibility could be that in fact there are 
strong correlations (ro • Yj) between flights along the cylindrical holes, so that the sum in (16.211) 
does not converge. 
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Figure 6.16: Integrated VACF R{t) for r = 1 .7, b = 0.0. 



0.545 1 
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Figure 6.17: Integrated VACF R{t) for r = 1.8, = 0.3. 
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Figure 6.18: Integrated VACF R{t) for r = 2.07, b = 0.6. Computed using an average over = 10^ initial 
conditions. 



Chapter 7 



Conclusions and future directions 



7.1. Conclusions 

We have studied geometrical, statistical and physical aspects of deterministic diffusion in three 
classes of bilUard models. 

We first discussed how best to estimate diffusion coefficients in billiard models from numerical 
data, presenting a method which estimates the width of the sampling distribution of the diffusion 
coefficient in terms of the rate of growth of the width of the distribution of the mean squared 
displacement. 

We apphed these numerical estimates in a 2D periodic Lorentz gas model to study the geometry- 
dependence of the diffusion coefficient, finding an unexpected quaUtative change in behaviour as 
one of the parameters is varied. We extended a previously-known random-walk approximation 
to our model, showing that there are two regimes where it can be applied. We then considered 
a Green-Kubo formula, the zeroth-order term of which is this random-walk approximation. We 
showed how to improve the justification of the derivation of this formula. We also made preUmi- 
nary investigations on the variation of the diffusion coefficient when the symmetry of the system 
is reduced. 

In the same 2D periodic Lorentz gas model we then studied the shape of position and distri- 
bution functions, which we know converge to normal distributions at long times. Using a method 
which we believe is more appropriate than the usual histograms in this case, we showed that at 
finite time the densities possess a fine structure which prevents them from converging pointwise 
to Gaussians. We showed how this fine structure could be understood in terms of the geometry of 
the billiard domain, giving an analytical expression for the fine structure function for both position 
and displacement distributions. 

Using these expressions, we showed that demodulating the densities by the fine structure gives 
functions which describe the two-dimensional density in a quasi-one-dimensional channel; in- 
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formation about these 2D density functions is otherwise difficult to obtain. This demodulation 
eUminates most of the very fine-scale oscillations, resulting in much smoother functions which 
seem to converge uniformly to Gaussian densities, strengthening the standard central Umit theo- 
rem. Nonetheless, in certain parameter regimes these underlying densities themselves possess a 
type of fine structure, now corresponding to a nearly-uniform distribution within traps. We also 
used the knowledge of the fine structure to give a physical picture of the rate of convergence in the 
central limit theorem. 

We considered the effect of imposing a non-constant distribution of particle speeds. We showed 
that for a general speed distribution with finite mean this does not affect the mean squared dis- 
placement. In the special case of a Gaussian distribution of velocities, however, we showed that 
the resulting limiting position distribution is skewed away from a Gaussian shape. Our analytical 
calculations of the resulting shape matched numerical results very well once the densities had been 
demodulated by the relevant fine structure function. 

We studied the extent to which similar statistical properties hold in a polygonal billiard channel 
for which there are very few rigorous results. Although chaotic and mixing properties are much 
weaker for this model, we showed that similar methods can be applied, giving evidence that normal 
diffusion can occur, in the sense that the mean squared displacement can grow asymptotically 
hnearly in time. We also confirmed that the central hmit theorem can be satisfied, but with a 
slower rate of convergence than for the Lorentz gas. We estabhshed in several particular cases that 
the existence of parallel scatterers in the structure results in anomalous diffusion, and conjectured 
that this is generally the case. We found a crossover from normal to anomalous diffusion when this 
geometrical configuration is approached. We were able to understand this qualitatively in terms 
of a continuous-time random walk model, although we found that the quantitative prediction of a 
simple version of that formahsm did not match numerical data well. 

We finally discussed to what extent results on two-dimensional periodic Lorentz gases can be 
extended to the more physically realistic three-dimensional case. We constructed a 3D periodic 
Lorentz gas with overlapping scatterers and showed that it has a finite horizon regime, in which 
it exhibits normal diffusion. We discussed how different types of holes in the structure affect the 
statistical properties, giving evidence that it may be possible to have normal diffusion even when 
corridors exist, provided they are small enough. 
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7.2. Future directions 

Our results point in several directions: 

• We hope that it is possible to find a better physical model of diffusion in the 2D Lorentz gas 
which can predict qualitative features of the geometry dependence of the diffusion coeffi- 
cient. 

• We would like to prove strong convergence of projected densities, as discussed in Sec. 14.11 
and App.lCl at least in a simple model such as the Arnold cat map ||Dor99ll . 

• It may be possible to develop a more quantitative version of the qualitative arguments given 
in Sec. 15.41 to model anomalous diffusion in polygonal billiards. In particular it would be 
interesting to derive (approximations to) the step length distribution in the continuous-time 
random walk model for the zigzag model, directly from the shape of the unit cell. 

• Our analysis of the fine structure in Chap. |4] may have implications for the escape rate for- 
malism for calculating transport coefficients (see e.g. IIGas98l ). where the diffusion equation 
with absorbing boundary conditions is used as a phenomenological model of the escape pro- 
cess from a finite length piece of a Lorentz gas; analyzing the fine structure in this situation 
could provide information about the validity of the use of the diffusion equation in that 
context. 

• Further investigation is needed of the effect of cylindrical holes in 3D models described in 
Sec. 16.61 It is important to establish if (and to what extent) it is possible for a 3D periodic 
Lorentz gas with non-overlapping scatterers to have normal diffusion. 

• The arguments discussed in Sec. l6.6l relating to the possibility of faster decay of correlations 
and free path distributions should extend to billiards in higher dimensions; in particular, 
we should in principle be able to obtain bounds on the decay of correlations and the free 
path distribution of hard-sphere fluids, by treating them as billiards in a high-dimensional 
phase space, as described in Sec. 11.2.31 An understanding of the shape of 'free' regions in 
the configuration space, together with arguments similar to those of Sec. 16.61 should give 
bounds on and/or estimates of the rate of decay of velocity autocorrelations and of the free 
path distribution, corresponding to the trapping effects referred to in ICYOOII . However, the 
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complicated geometry of the high-dimensional phase space IIChe97ll means that this would 
be difficult to implement. 

We expect that bilhard models will remain of interest to mathematicians and to physicists in 
the future. 



Appendix A 



Convergence of rescaled solutions of the diffusion 
equation to a Gaussian 



We show that the rescaled solution of the diffusion equation starting from an initial density which 
decays sufficiently rapidly at infinity converges to a Gaussian shape as f ^ oo. 

A.l. Pointwise convergence 

Let pt be the solution of the diffusion equation with initial condition po which is a density, i.e. 
which satisfies po ^ and / po (y) dy = I. We also assume that po : M ^ M+ is a piecewise 
continuous function and decays sufficiently fast at infinity. By translating the coordinate origin if 
necessary, we further assume that the centre of mass of the initial condition is fixed at the origin: 
/3'Po(3')d3' = 0. 

Define the rescaled solution pt{x) := ^/tpt{x^/t). Then 



For fixed x, the argument of the exponential in the integrand in (IA.2I ) tends to as f — > oo. 
Hence the integrand tends to pQ{y) and it is bounded above provided po{y) decays sufficiently fast 
at infinity, for example if it is exponentially locaUsed in the sense that 




(A.l) 




(A.2) 



Po{y) ^ e 



■K\y\ 



(A.3) 



for some constant K>0. 
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The Lebesgue dominated convergence theorenj^ IIRud76ll then implies that for fixed x, 



V471D 



,-x^/4D 



(A.4) 



i.e. the rescaled solution of the diffusion equation tends pointwise to a Gaussian with variance 2D. 
In fact the convergence is uniform, as follows. 



A.2. Uniform convergence of rescaled density functions 

Consider the difference 

\Pt{x) - g2D{x)\ , 

where gfji is the Gaussian density with mean and variance given by 
The difference is then given by 



\ptix) - g2D{x) 



-x^/AD 



po{y) exp 



xy y'- 



lD^ft 4Dt 



dy-\ 



We may bring the —1 term inside the integral, since / po{y)dy = 1 by assumption. 
We expand the exponential in a Taylor series for t large: 



exp 



xy y 



2 1 



2Dy/t 4Dt 



n=i 

xy 



xy y- 



2Dy/t 4Dt 

2 2 2 

y X y 



xy-' 



2D^t 4Dt ' 8D2f 8D2f3/2 ' 32D^t^ ^^^3 



+ ■ 



Assuming that it is permissible to integrate term-by-term, we have 

\pt{x) - g2D{x)\ = g2D{x) ^^JyPo{y)dy+(^-^^-^^ //po(y)dy 



i „-x^/4d{x^)0 ( X 



t \ 8D2 4D 



(A.5) 



(A.6) 



(A.7) 



(A.8) 
(A.9) 



(A. 10) 
(A.ll) 



^If (/„) is a sequence of measurable functions such that f„(x) f{x) almost everywhere as n ^<x> and there exists 
an integrable function g such that |/n(j^)| g(x) almost everywhere, then lim,j^oo/ fn = J f- 
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since / ypo(y)dy = by our choice of coordinates and (x^)o = / y^po{y) dy by definition. 

Since the decay of e^''^/^^ is faster than the growth of any polynomial in x, the terms in x in the 
above are bounded. If also the term written as 0{t^^^^) is bounded, then we have the following 
estimate for large t which is uniform in x (i.e. independent of x): 

\Pt{x)-g2D{x)\^j, (A. 12) 

for some constant C. Thus ||p, — ^2Z)|L ^ C/t for sufficiently large t, where the uniform norm is 
defined by 

||/L:=sup|/(x)|. (A.13) 
Hence the rate of convergence of the rescaled density to the limiting Gaussian is C'(?^'). 

A.3. Rigorous proof of uniform convergence 

We now give a rigorous proof of the result on uniform convergence for which a heuristic argument 
was given in the previous section. We adapt a method from MMilOOl . where convergence to a 
Gaussian was considered without rescaling; see also IIMBOBI and references therein. We remark 
that a faster rate of convergence was obtained in MMilOOl by choosing a different time origin for 
the Green function, but this does not work in our situation due to the rescaling. We set D = 1 for 
simplicity (e.g. by rescaling time). 

Theorem A.3.1. Let the density po be piecewise continuous and such that the first two mo- 
ments J ypo{y)dy and f y'^ Po{y)dy exist. Let pt be the solution of the diffusion equation with 
diffusion coefficient D = 1, starting from the initial condition po. Then the rescaled density 
Pt{x) := \/t Pt{x\/1) converges uniformly to the limiting Gaussian G'{x) := -^^z^^l'^ as t ^ oo, 
with rate of convergence 0{t^^). 

Proof. Let the initial condition be h{x) := p{0,x). Then the solution at time t is given by the 
convolution Pf = h*G', so that taking Fourier transforms gives Pt{k) = h{k)e^'^^', where the 
second term is the Fourier transform of the Green function G'{x). 

Let (x) : = Pt (x) — G' (x) be the error of the solution at time t compared to the Green function, 
and let 

E'{x) := ptix) - &{x) = VtE'{xVt) (A.14) 
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be the error of the rescaled solution from the Umiting Gaussian. 

Then the Fourier transform of the error is 

E\k) = {h{k) - (A.15) 

By Taylor's theorem with remainder we can expand h{k) as 

h{k) = h{0) + kh'{0) + ^k^h"{c), (A.16) 

for some c = c{k) G [0,A:]. But ^(0) = 1 and ^'(0) = — i(x)o = 0, by our choice of coordinates. 
Furthermore, 

oo 

h"{c) = J {-ixf h{x)dx, (A.17) 

— oo 

so that 

oo 

\h"{c)\f^ J x^h{x)dx<oo, (A.18) 

— oo 

since we assumed that the second moment of h exists. Thus 

^CA:V*'', (A.19) 

for some constant C. 



We now need to convert to the rescaled functions. We have 

oo 

E'{k) = j e-''^VtE'{xVt)dx (A.20) 

— oo 

= E'{k/Vt) (A.21) 
after making the change of variables y = x\ft. 

Hence 

n\k)\ = \E\kl^t)\^^k\-^ . (A.22) 
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Reverting to real space using the inverse Fourier transform, we have 



\E'{x)\ 



2k 



E'{k) 



dk ^ 



Int 



k^e^'"' dk : 



C 

2Kt 



(A-23) 



where C is a constant independent of x. Thus pt converges uniformly to as f — > oo, and the size 
of the error is 



\E'\ 



\p'-&\\=0{t 



(A.24) 
□ 



The numerical results reported in Sec. 1 1.3 1 provide evidence that this upper bound t ' for the 
asymptotic rate is in fact the actual decay rate. 

A.4. Convergence of distribution functions 

We have shown that the convergence of the rescaled density functions is uniform in the case of 
the diffusion equation. However, for the diffusive dynamical systems considered in this thesis, 
rescaled density functions do not usually converge even pointwise to a limiting Gaussian distribu- 
tion; rather, we must consider convergence of the (cumulative) distribution functions. 

Let Ft be the distribution function at time t, given by 



Ftix) := I p,(x')dx', 



(A.25) 



and N be the distribution function of the limiting Gaussian distribution, so that 



N{x) := 



(A.26) 



Then Pt{x) = F/{x) (where the prime denotes differentiation), so that Pt{k) = ikFt{k). Follow- 
ing the same type of argument as for the density functions gives 



\\Ft-N\\^^^^=0{t 
for sufficiently large t and some constant C' . 



-1/2N 



(A.27) 
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The rate t^^l^ is known to be the fastest rate in the central limit theorem for independent and 
identically distributed random variables IIFelTll . and is also the maximum rate that we find in 
billiard models: see Sec. 14.61 
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Suspension flows 



Since rigorous results on the billiard flow are usually proved by using the fact that it is a suspen- 
sion over the bilhard map T, under the free path function T, we recall the definition of suspension 
flows and some key properties. A clear recent reference on limit theorems for suspension flows, 
with strengthened versions of relevant theorems, is IIMT04II . 

B.l. Definition of suspension flows 

We follow closely the definitions in Cornfeld et al. IICFS821 Chap. 11]. 

Let {X,B, v) be a measure space with an automorphism T , and let r : X — > M_|_ be a function 
such that /x?'dv < °°. 

Define the space under the roof function r by 

Y ■=X'' ■={{x,t):xeX, Q^t<r{x)}. (B.l) 

We assign a sigma-algebra B onY hy taking as measurable sets the measurable subsets of X x M 
which belong to Y , and we put 

^:=lvx£, (B.2) 
r 

where 

r ■= Ev [r] := J r{x)dv{x) (B.3) 

X 

and i is Lebesgue measure on M+, so that 

jU(A) = jJJ dv{x)dt. (B.4) 

A 

This gives a measure space {Y,B,v) with the normalisation IJ.{Y) = 1. 

We visualise the space 7 as a subset of the Cartesian product X x M+, as in Fig. IB. II We then 
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define the flow under the roof function r by flowing vertically from (x, 0) at unit speed for a time 
r{x), until we hit the roof function at (x,r(x)), when we instantaneously jump to (r(x),0). This 
corresponds to identifying the points (x,r(x)) and (r(x),0), in which case we can write the flow 
as 



computed using the identification: see IICFS821 [Gas98l for explicit expressions. 

The following result shows that any flow satisfying certain conditions can be viewed as a 
suspension flow, which is technically simpler to study. For a proof, see IICFS821 pp. 295ff]. 

Theorem B.1.1. Any flow V without flxed points on a Lebesgue space {M,B,n) is measure- 
theoretically isomorphic to a suspension flow (also called a special flowj. 

B.2. Ergodicity of suspension flows 

The following theorem is stated in IIChe02[ Sec. 4], but I could not find a proof in the literature so 
I give one here 3- 

Theorem B.2.1. Let be a suspended flow over the transformation T : X ^ X, under the roof 
function r : X ^ Then V' is ergodic if and only ifT is ergodic. 

Proof. Firstly suppose that the map T is not ergodic. Then there exists an invariant set A for the 
map which has non-trivial measure, i.e. /i(A) 7^ and / 1. Define the set B by 



see Fig. lB.ll Then B is an invariant set for the flow with measure v{B) which is non-trivial, i.e. 
v(B) / and v(S) / 1. Hence the flow V' is not ergodic. 

Conversely, suppose that the flow is not ergodic, so that there is an invariant set B for the flow 
with non-trivial measure. Then the set B must be of the form (IB. 61 ). since otherwise B would not 
be invariant. Hence the set A defined by projecting B down to X has non-trivial measure, so that 
the map is not ergodic. □ 

Thus if a flow can be expressed as a suspension over some map, then the ergodicity of either 
implies the ergodicity of the other. We need the following application of this in Sec. 13.5.31 

' I would like to thank Peter Walters for showing me the idea of this proof. 



V'{x,s) = {x,s + t) 



(B.5) 



B := {(x,0 : X G A, ^ f < r(x)} = |J{x} x [0,r(x)); 



(B.6) 



B.3. Central limit theorem for suspension flows 
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Figure B.l: Invariant sets A for the base transformation and B for the suspended flow. A does not need to 
be connected in general. 

Corollary 6.2.2. The torus-boundary map, which maps one intersection with the torus boundary 
into the next, is ergodic. 

Proof. We consider the bilUard dynamics on the torus. Sinai IISin70ll proved that the billiard map 
T, which takes one collision with a scatterer to the next, is ergodic. Hence the billiard flow <I> 
is ergodic, by Theorem IB. 2.11 since it is a suspension over the billiard map under the free path 
function T which gives the time taken from one collision to the next. 

However, the billiard flow can also be viewed as a suspension over the torus-boundary map, 
under the trap residence time function p . Using Theorem IB.2.11 again, we find that the torus- 
boundary map is ergodic. □ 

B.3. Central limit theorem for suspension flows 

The standard method to prove statistical properties for suspension flows is via the above construc- 
tion, relating them to statistical properties of the Poincare map. 

Let fr '■ = lo f ° df. We say that the central limit theorem (CLT) is satisfied for fx if 




(B.7) 



for some normal random variable Z. 
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B. Suspension flows 



Define F{x) := J^^' f{x,u)du. The following theorem is proved in IIMT04I1 under certain 
technical conditions. 

Theorem B.3.1. Suppose that F and r both satisfy the CLT. Then f satisfies the CLT. If the CLT 

for F has variance of ^ 0, then the CLT for f has variance = of /f, where f := Im^^^ 
is the mean of the roof function. 



Appendix C 



Convergence of projected densities 



C.l. Densities and the Perron-Frobenius operator 

We consider a flow ^ : Ai —>^ M on the phase space M. Let /Iq be a measure describing the 
distribution of initial conditions at time t = 0. After time t, this has evolved to the pushed-forward 
measure jXt defined in Chap. |2] 

We say that /Iq is absolutely continuous with respect to /i, denoted Hq <^ n, if IJ,{A) = ^ 
Ho{A) = 0, i.e. if Hq is not concentrated on any sets of zero /i -measure. In this case, the Radon- 
Nikodym theorem |Roy88| shows that there exists a unique non-negative function /o G L^(A^) 
such that 

Ato(A) := Jfodn := J fo{x)dn{x) for all A G B. (C.l) 

A A 

We call /o the density of /Iq with respect to /x. 

If /X is measure-preserving, so that ;U,(A) = /Xo(<I>^'(A)) for all A G i3, then takes sets of 
measure to sets of measure 0. Hence if Hq <^ jj. then also /i, <C /X, so that Ht also has a density, 
which we denote by /(. The map P' given by 

which describes the time evolution of phase space densities is called the Perron-Frobenius oper- 
ator IILM941|KH95II . In the case of invertible, measure-preserving transformations, we can write 
an expUcit formula for the time-evolved density IIKH951 Chap. 5] : 

f,{x)=h{^-'{x)). (C.3) 

This shows that in this case the density just gets 'moved around' ; nonetheless if the flow is mixing 
then this 'moving around' occurs in such a way that the density gets spread out over phase space, 
as follows. 
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C. Convergence of projected densities 



C.2. Mixing and weak convergence of densities 



Recall that the flow 



is mixing with respect to the invariant measure if 



At (A) At (B). 



(C.4) 



We can re-express this in terms of functions as follows: 



/ 



llo'(A)lBdAt ^ / UAdAt / UfidAt 




(C.5) 



M 



M M 



By an approximation argument the following theorem relating mixing to convergence of den- 
sity functions can then be proved: see e.g. IILM941 p. 73]. 

Theorem C.2.1. Let {M,B,lJ.) be a probability space (i.e. such that }x{M) = I), : M ^ M 
a measure-preserving flow, and P' the Perron— Frobenius operator corresponding to <!>'. Then 
is mixing if and only if (P'f) is weakly convergent [see below] to Ij^ for all f D, where 
D:={feL\M):f^O}andlM (x) = 1 for all xeM. 

(Note that / ^ does not strictly make sense for an function, whose values can be changed 
on a set of measure without affecting the function. The meaning is that it is possible to find a 
representative of the equivalence class for which / ^ 0, or equivalently that f{x) ^ for almost 
allx£M.) 

We now define the notion of convergence appearing in the theorem. 

Definition C.2.2. A sequence (/„)„gn» /« ^ L^, converges weakly to f £ L^, denoted /„ ^ /, if 
and only if 



lim(/„,g) = (/,^) for all LP', 



(C.6) 



where 




(C.7) 



M 



Here, is the dual space of with 



— I = 1 forl</7^oo 

p p' 



(C.8) 
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and p' = oo when p = I, where L°° is the space of essentially bounded functions, i.e. functions 
which are bounded except on a set of measure 0. 

Corollary C.2.3. is mixing if and only if 

{P'f,g) ^ (1^^,^) = Jgdu ast^ oo, (C.9) 
M 

for all bounded functions g : M ^M. 

C.3. Weak convergence of projected densities 

In general we cannot have stronger than weak convergence of densities in phase space to the 
invariant density. For example, it is proved in IIGLS98II that there cannot be convergence in the 
norm for any q ^ I, since the measure of the set with density in any interval is conserved. 

However, we might expect that marginal densities obtained by projecting onto lower-dimensional 
subspaces of the phase space may be able to converge more strongly ||Dor99l| . Here we prove that 
they converge weakly; in the next section we discuss the question of strong convergence. 

Consider for concreteness the 2D periodic Lorentz gas, with coordinates {x,y,d) in the phase 
space M = QxS^. Consider an initial distribution in phase space given by the density /o : — > R 
with respect to normalised Liouville measure d^u := j^^^dxdydd, where dx is the differential of 
Lebesgue measure in the A-direction. This density evolves in phase space via the Perron-Frobenius 
operator P', defined as above. 

We define the measure /i' by d/i' := y^dxd)', i.e. normalised Lebesgue measure on Q, and 
projected densities 0, : 2 ^ M on the configuration space Q by 

271 

(j)r:= J{P%)dd= J f{x,y,d)dd, (C.IO) 
5' 6=0 

setting ft := P'fo- Then 0, G ^^(m') and 0, ^ 0, so that <pt is a density. 

We wish to show that the (0,) converge weakly to 1q, the constant invariant density on Q, with 
respect to the measure /i'. We distinguish when necessary the measure over which we integrate 
by writing 

{f,8)n ■■= J fgd^l (C.ll) 
M 
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C. Convergence of projected densities 



for the inner product of / and g with respect to the measure jj., although in principle the measure 
is implicit in the domain of the functions / and g. 



Let 7 G L°°{Q) be a bounded function on Q. We want to show that 



{^t,r)ti'^{'^Q,r)fi' = Jr^^' as? 



(C.12) 



by relating the left hand side to objects in phase space. 



We have 



1 

"W\ 



2n 



ft{x,y,d)de 



Yix,y)dxdy 



1 



1 el 



ftix,y,d)g{x,y,e)dxdydd, 



QxS' 



where we define g by 

g{x,y,e) :=Yix,y)h^{e), 
so that g is constant on fibres over Q. Thus 

QxS' Q 

But yGL°°{Q) was arbitrary. Hence 0^ ^ Ig, as required. 



(C.13) 

(C.14) 
(C.15) 



(C.16) 



(C.17) 
(C.18) 



C.4. Convergence of ID distributions in billiards 

The above constitutes a general method for such proofs. We now consider the special case of 
projecting down from densities on Q to densities in one coordinate direction, as in Chap.lH 



We define the measure v' on the x-space to be Lebesgue measure. We then define ID projected 



C.4. Convergence of ID distributions in billiards 
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densities in the x-direction by 



i 

¥t{x) := J (l>tix,y)dy. 

y=0 



(C.19) 



(We could write Yt '■= /si <l>t dy, since in fact the x and y coordinates are defined on a torus.) 



To study weak convergence of the xi/t, consider an arbitrary bounded (i.e. L°°) function p. To 
mimic the previous proof, we wish to define y. Q^M. such that 



But 



{Wt,p)v = {<l>t,y)^i'- 



1 1 

{Wt,p)v= J J (l>tix,y)dy 



x=0 \y=0 



p{x)dx 



JJ (l>t{x,y)Y{x,y)dxdy 



x,y 



if we set 

r{x,y) :=p{x)lH(^^{y). 
The subtlety here is that different fibres have a different amount of associated measure. 



(C.20) 



(C.21) 
(C.22) 

(C.23) 



With the above definition of 7, we have 



{¥t,p)v' = {<l>t,7)n' 

^ {h^r)^l' = J rdiU' = / P{x) J ]l//(x)(j)dj 

Q X 
= J p{x)h{x)dx= {h,p)v'- 



dx 



Thus 



(C.24) 
(C.25) 

(C.26) 
(C.27) 



with respect to the measure v'. 
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C. Convergence of projected densities 



A different point of view is to consider the canonical projection 

k:Q^S^; {x,y)^x. (C.28) 

Denoting by v tiie push-forward of the measure fJ.' under this projection, we have 

v(A) := [7r*(Ai')](A) -^'{K-'iA)) = J ]l,eA]l(.,,)eedM' (C.29) 

Q 

= J l{;cGA}U{(.v,v)ee}dM' (C.30) 
Q 

= J J dydx= J h{x)dx. (C.31) 

xeAyeH{x) xeA 

Thus the natural measure v on the x-space has density h with respect to Lebesgue measure, so that 
we could equally look at the density v/'(x) := i/a(x)//j(x) with respect to the measure V and say 
that 

(v/,p)v^(l,p)v, (C.32) 
i.e. that xj/ converges weakly to the constant density with respect to the geometrical measure v. 

C.5. Stronger convergence of projected densities? 

We would like to prove that projected densities converge strongly, e.g. in or even uniformly. 
This cannot be true in general, since if we project along the stable direction then we do not ob- 
tain any smoothing effect: see the discussion of the baker and cat maps in ||Dor991 . We expect, 
however, that if we avoid this special direction then we should get strong convergence. We are not 
aware of any rigorous results on this, even for relatively well-understood systems such as the cat 
map. However, Sinai proves in IISin94[ Chap. 18] that densities projected to unstable manifolds 
converge pointwise; see also IIGLS98I . 



Appendix D 



Weak convergence of measures in path space 



D.l. Measures on path space 

We recall the definition of convergence in distribution of the random path % to the Wiener process, 
a convergence of measures on the space of paths. 

Define the accumulation function by 

t 

Sr{-):= lfo^\-)ds, (D.l) 


so that St : M ^ R. We denote elements of M by ft). 
Define a rescaled process Wj by MCYOOl 

Wri..co):=»^^, (D.2) 

where T > 0, t = sT , and the path is parametrised by 5 G [0, 1]. 

For fixed T, Wt{-',(o) : [0, 1] ^ is a continuous path in M^. In the case of diffusion in 
billiards, we have 5,(ft)) = jQ\{s;(o)ds = Xf(ft)) — xo{(o), where ft) = {q,v) denotes the initial 
condition. 

The continuous random function Wt induces a probability measure Pj on the space C([0, 1]) 
of continuous paths from [0, 1] to Ai, via 

Pt{A):=h{(0£M:Wt{-;(o)£A), (D.3) 

where A C C([0, 1]) is a subset of the space of continuous paths which is measurable with respect 
to the Borel (7-algebra on the metric space C([0, 1]) with metric 

d{f,g):= sup \f{x)-g{x)\. (D.4) 

A-e[0,l] 
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D. Weak convergence of measures in path space 



D.2. Weak convergence of measures on metric spaces 

The notion of convergence in tiie standard central limit is convergence in distribution of the 
rescaled distributions to a normal distribution, which can be expressed in terms of pointwise (and 
uniform) convergence of the rescaled distribution functions to a normal distribution function on 
R" IIB 11681 . This definition in terms of distribution functions cannot be directly generalised to 
convergence in path spaces, which we require here, but it is equivalent to the following notion of 
weak convergenci^ of measures, which does generalise to arbitrary metric spaces MB 116811 . 

Let {X,B) be a metric space together with the a-algebra of Borel sets on it. Let {Pn)neN and P 
be probability measures on {X,B). Then f„ converges weakly to P, written P„ P, if and only if 

I fdPn ^ J fdP, (D.5) 

X X 

for all bounded, continuous functions / : X ^ M ||Bil68i . There is an equivalent formulation in 
terms of sets, reflecting the duality between measures thought of as the dual space of functions 
and measures as set functions. Namely, P„ P if and only if 

Pn{A) P{A) for all A such that P{dA) = 0, (D.6) 

where dA:=A \A° is the boundary of the set A, i.e. the set of points which are limits of sequences 
of points in A and limits of sequences of points outside A. (Since dA is a closed set, it belongs to 
the Borel a-algebra B, so that P{dA) is defined.) 

Sufficient conditions for this weak convergence are IIB 116811 : (i) the finite-dimensional distri- 
butions converge; and (ii) this convergence is tight. Property (i), in the case of convergence to 
Brownian motion, is a multidimensional central limit theorem IIChe95l . In Sec. 12.4.2] we reformu- 
late it using the more intuitive notation IIBle92ll 

,,[0,,| (D.7) 

V' 

for the rescaled process, where %{s) =Wt{s,-) : M. — > M. Property (ii) means that for all e > 0, 
there is a compact set^ = ^(£) such that Py (A') > 1 — £ for all T; this prevents mass from escaping 
to infinity ||Bil68l . 



'Weak convergence' in probability theory is close to weak-* convergence in analysis IRWOOl Section II.6]. 
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